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Abstract 


We study the short pulse dynamics in the deterministic and stochastic environment 
in this thesis. The integrable short pulse equation is a modeling equation for ultra- 
short pulse propagation in the infrared range in the optical hbers. We investigate the 
numerical proof for the exact solitary solution of the short pulse equation. Moreover, 
we demonstate that the short pulse solitons approximate the solution of the Maxwell 
equation numerically. Our numerical experiments prove the particle-like behavior of 
the short pulse solitons. Furthermore, we derive a short pulse equation in the higher 
order. 

A stochastic counterpart of the short pulse equation is also derived through the use 
of the multiple scale expansion method for more realistic situations where stochastic 
perturbations in the dispersion are present. We numerically show that the short 
pulse solitary waves persist even in the presence of the randomness. The numerical 
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schemes developed demonstrate that the statistics of the coarse-graining noise of the 
short pulse equation over the slow scale, and the microscopic noise of the nonlinear 
wave equation over the fast scale, agree to fairly good accuracy. 
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Chapter 1 
Introduction 

1.1 Introduction 

Physical phenomena whose states change over time are modeled by different mathe¬ 
matical techniques. Often, one is interested in the behavior of the physical system in 
time, or hnds it necessary to study of the dynamical behavior of the physical system. 
Modeling equations are generally given in partial-differential equation (PDE) forms 
[1]. The wave equation, the heat equation and the Schrodinger equation are only a 
few famous partial differential equations that we encounter in physics and the applied 
sciences. If the modeling equation for a given dynamical system is integrable, such a 
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system is known as an integrable system. These types of models are not only used in 
the natural sciences and engineering, but also in social sciences such as hnancial and 
economic forecasting and environmental modeling [2]. Not all physical phenomena 
can be modeled in a deterministic manner. Deterministic dynamical systems need 
not be linear. Modeling equations may also appear in nonlinear differential forms [3] 
whose solution may not be found analytically. Numerical techniques are then utilized 
to study such dynamical systems mi- On the other hand, some other phenomena 
appear to be highly stochastic. Stochasticity may arise from an internal or external 
mechanism in the system and may have some visible impact in the evolution of the 
system state. Under the influence of randomness, the system is not deterministic any¬ 
more and it requires a probabilistic and statistical approach to study the evolution 
of the system [0]. Governing equations for such systems can be written in stochastic 
differential form [7] or in corresponding Fokker-Planck equation form [8]. 

Interaction of light with matter is an important physical phenomenon that has 
played a central role in the recent history of technological advances. Light is an elec¬ 
tromagnetic wave. A wave is simply a quantity that varies with time and position. It 
is a periodic disturbance with many oscillations, and propagates through a medium 
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such as empty space or an optical fiber [9]. A wave is generally modeled mathemat¬ 
ically by using partial differential equations (PDEs) with a wave value an 

independent variable time t and one or more independent spatial variables x. Waves 
can fall under two categories: linear waves and nonlinear waves. Linear waves are 
those modeled by linear equations. The classical wave equation in three dimensions 
is very well-known example: 


1 d‘^u 

dU 


Vu, 


( 1 . 1 ) 


whose solution is a linear wave. Linear PDEs are easy to solve as opposed to nonlinear 
PDEs. Since the superposition principle applies to linear waves, linear combinations 
of the simple solutions can be used to form more complex solutions for complicated 
linear problems. However, the principle of superposition does not apply to the sec¬ 
ond category of waves, namely, nonlinear waves. Nonlinear waves are described by 
nonlinear equations. The nonlinear Schrodinger equation (NLSE), the Korteweg-de 
Vries equation (KdV) and the sine-Gordon equation (sG) are some of the very well- 
known and well-studied examples of integrable nonlinear equations whose solutions 
are solitary waves, also called solitons. The inverse scattering theorem is used to 
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solve integrable nonlinear PDEs Although the nonlinear KdV and sG equa¬ 

tions have exact solitary wave solutions, there are some other types of nonlinear wave 
equations whose solutions display singularities or those decaying in time and space 
(smooth dispersive solutions). 

The discovery of optical fibers is a milestone in modern day communication tech¬ 
nology and has been used in many optical systems since then. Compared with the 
electrical transmission system based on copper wires, fiber optic communication uti¬ 
lizes huge amount of data to transmit signals. That serves as the main reason why 
optical fibers based communication systems have played a major role in the advent 
of the information age. The set of Maxwell’s equations is an elegant description of 
light and matter interaction, and can be reduced to a single nonlinear wave equation 
for studying optical phenomena. However, the exact solution of the nonlinear wave 
equation is unknown as of yet. Nevertheless, the wave equation can be approximated 
to simpler forms expressed as nonlinear partial differential equations to describe light 
propagation in optical fibers. 

The asymptotic multiple scale expansion of the nonlinear wave equation yields 
the cubic nonlinear Schrodinger equation m- The assumption that is made in the 
derivation of the NLSE is that the pulse width is large in comparison to the carrier 
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wavelength [12]. The NLSE is the governing eqnation for pulse propagation in opti¬ 
cal hbers and possesses exact solitary wave solutions 0. As the name suggests, the 
NLSE is just the nonlinear version of the famous quantum mechanical Schrodinger 
equation. It can also be applied to some other physical phenomena such as hydrody¬ 
namics and quantum condensates m. 

The short solitary waves, or short pulses, are used to study the nonlinear effects 
in optical hbers due to their unique solitonic properties. The width of those short 
pulses generally ranges from 10 nano-seconds to 10 femto-seconds [I5]. As the pulse 
shortens even further, the NLSE fails to be a good modeling equation [T21 dSl [Il|. 
A nonlinear partial differential equation derived by T.Schafer and C.E. Wayne m 
may be used to model the propagation of ultra-short pulses [T^ EH] in nonlinear cubic 
media (optical hbers) and is given, up to a scale transformation, as 


Mri = M + . 


( 1 . 2 ) 


This equation is called the short pulse equation (SPE). The SPE is integrable [TU] . 
and possesses exact one-soliton I2D1 and multi-soliton es solutions. There has been 
an intensive research on the SPE over the past few years due to its integrability and its 
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exact soliton solutions. The Hamiltonian structures [22l |23] and conserved quantities 
PI of the SPE are found. The periodic solutions are also obtained pg. The vector 
short pulse equation (VSPE) is introduced [26] and the integrability of the VSPE is 
studied as well [27|. The regularized short pulse equation (RSPE) is derived [2H|- The 
existence of multi-pulses [29] , traveling waves, and solitary wave solutions [2HI EO] of 
the RSPE are studied. 

The study of the short pulse equation in both deterministic and stochastic media 
is the main focus for the rest of the chapters. We present Maxwell’s equations and the 
nonlinear wave equation in lieu of the optical phenomena discussed in chapter two. 
The discussion of the multiple scale expansion method and the NLSE is presented in 
chapter three. The derivation of the SPE followed by a numerical analysis is provided 
in chapter four. We also derive a higher order SPE in this chapter. In chapter Eve, we 
investigate SPE solitons and their particle-like properties. The following chapter deals 
with a relatively new area concerning stochasticity. The derivation of the stochastic 
SPE, and the discussion on coarse-graining noise are given in chapter six as well. 
Since the numerical analysis is an important portion in our work, we reserve an entire 
chapter for numerical methods employed in our numerical experiments. The closing 
chapter provides prospect for research pursuits regarding short pulse dynamics. 


Chapter 2 


Pulse Propagation in Nonlinear 
Cubic Media 


The theory of electromagnetic wave propagation in dispersive nonlinear media has 
played a major role in the advent of the 20th century’s information technology. In 
this chapter, we will derive a basic equation that governs propagation of optical pulses 
in nonlinear cubic media such as single-mode hbers in one dimension. This equation 
is the starting point of our discussion in this thesis. 
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2.1 Maxwell’s Equations 

The theory of wave motion is an important mathematical model in many areas of 
physics and engineering. A large number of real world applications can be explained 
using the solutions of the wave equation. We look at the wave theory in the perspective 
of light and obtain a one dimensional model associated with the light propagation in 
nonlinear cubic media. Mathematically, the basis of wave theory is the wave equation 
which is a second-order partial differential equation. One can derive a wave equation 
either in a linear form or nonlinear form depending upon the polarization of the 
material [an [12]. The linear wave equation is no longer sufficient to describe the 
propagation of light in a medium when the intensity of light becomes large enough. 
In such situations, light waves interact with one another and with the optical medium 
leading to the emergence of nonlinear effects. These nonlinear phenomena require an 
extension of the linear theory, and a nonlinear response of optical materials to the 
optical excitations must be taken into account. In the one-dimensional case, the 
displacement of the wave is assumed to be a scalar function and the pulse is a scalar 
wave. One must refer to Maxwell’s equations to grasp a deep understanding of wave 
motion in a variety of situations. Maxwell’s equations are the governing equations 
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for the propagation of optical pulses in nonlinear media and are given as a set of 
equations in differential form in three dimensions as [32] 


V X E 

V X H 


^ dB 


VD= p 
V-B= 0 


( 2 . 1 ) 


where p is charge density, J is current density. The magnetizing held or magnetic 
held H and electric held E are related with the corresponding magnetic hux density 
or magnetic induction B and the electric displacement held D through 


B = /xoH + M 

(2.2) 

D = cqE + P 

(2.3) 


where P is electric polarization, M is magnetic polarization, eo is the electric constant 
and Pq is the magnetic constant. Let’s briehy mention the physical meanings of 
these quantities before deriving the wave equation from Maxwell’s equations. The 
electric constant Cq is also called the permittivity of free space and describes the 
interaction of free space with an electric held. Therefore, permittivity is a physical 
quantity and relates the electric charge to the mechanical quantities such as force. 
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Similarly, the magnetic constant /xo is also a physical constant which is sometimes 
called the permeability of free space. It is the degree of magnetization that vacnnm 
can obtain in response to an applied magnetic held. As for the electric polarization, it 
is simply the average electric dipole moment per nnit volnme and can be formnlated as 
P = Po where pi is the dipole moment of the ith molecnle and N is the nnmber 
of molecnles per nnit volnme. Notice that the average of the snm of the dipole 
moments of all molecnles comes from the fact that N is the statistically averaged 
large nnmber of molecnles. In general, the polarization vector P depends npon the 
local valne of the electric held strength E nonlinearly. The total polarization indnced 
by electric dipoles in a material as a response to an applied electric held satisties the 
general form [31] and can be written as a power series in the electric held strength 

P = eoix^^lE + . e2 + ;^(3):e3 + ...) (2.4) 

where is the jth order susceptibility. The nonlinear susceptibility is a quantity 
that is used to determine the nonlinear polarization of a medium and is a very useful 
quantity when describing nonlinear optical phenomena. In a similar fashion, M is 
the magnetic polarization and is simply an average magnetic dipole moment per unit 
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volume. There is no further explaination required for this concept because it is zero 
for nonmagnetic materials. By manipulating Maxwell’s equations, we can obtain 
a number of the properties of light such as the relationship between the E and B 
helds. Maxwell’s equations require light to be a transverse wave, i.e., the vector 
displacements E and B are perpendicular to the direction of propagation k. It will 
be seen in the next section that we choose E and B helds based on this condition. 
In light of the physical meanings of these concepts. Maxwell’s eqnations can now 
be rearranged to display explicitly the time and coordinate dependence of the wave 
amplitude. 

2.2 The Nonlinear Wave Equation 

To derive a wave equation that describes pnlse propagation in nonlinear cnbic media 
snch as optical hbers in one dimension, we choose the electric held E in the direction 
and the H held in the y direction snch that 


E = u{x, t)k, H = H{x, t)j 


(2.5) 
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where u{x,t) and H{x,t) are the magnitudes of electric field and magnetic fields re¬ 
spectively. Remember the choice of perpendicular electric and magnetic fields to each 
other and to the direction of propagation is the requirement of Maxwell’s equations. 
It should be noted that p = 0 and J = 0 in the absence of free charges, and M = 0 
for nonmagnetic media. Because we consider a nonmagnetic nonlinear cubic medium 


with no free charge in our model, we set p = 0 and J = 0, and M = 0 in (2.1) 


and (2.2) respectively. Once we substitute (2.5) in (2.1), we reduce the set of three 


dimensional equations (2.1) to a set in one dimension 


du 


dH 


dx dt 
dH d{eoU + p) 
dx dt 


( 2 . 6 ) 


where p = p{x, t) is the magnitude of polarization along the 2 ;-direction. To obtain an 
equation for the description of pulse dynamics in terms of the magnitude of the applied 


electric field, we take the derivative of the first equation in (2.6) with respect to x 


and the derivative of the second equation in (2.6) with respect to t. Subsequently, 


elimination of the H terms and a basic manipulation in the set lead to the wave 
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equation 


d'^u 1 d'^u 


d'^p 


dx'^ dU dP 


(2.7) 


where c is the speed of light in vacuum and the relation c = 1/ ^/JM^ is used. Equa¬ 


tion (2.7) is the standard wave equation that describes the propagation of linearly 


polarized light in one dimensional medium with u = u{x,t) being the magnitude of 
the applied electric field and p being the polarization of the medium in response to 


the electric field. Equation (2.7) is not a good representation of wave motion because 


it contains both magnitudes of the electric held and polarization. However, we al¬ 
ready know from the previous section that polarization can be expanded in terms of 
the electric held strength. Such a relation between polarization P and electric held 
E can be used to obtain an equation in terms only of the magnitude of the electric 
held. The total polarization induced by electric dipoles in a material as a response 


to an applied electric held satishes the general form (2.4). The main contribution to 


polarization comes from the linear term. For media that have inversion symmetry, the 
quadratic term in the expansion of the polarization vanishes. Nonlinear cubic media 
(for example, silica hbers) do exhibit such symmetry and no quadratic susceptibility 
contributes to nonlinear ehects. Hence, the third order susceptibility is the origin of 
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the lowest-order nonlinear effects in optical fibers. Nonlinearity is very important for 
the application of optical data processing, and the physics of nonlinear effects can be 
extracted by studying the response of the applied optical field to the medium. Since 
the contribution of higher order nonlinearities to the total polarization is negligibly 
small, we truncate the series at the third order of susceptibility, and therefore consider 
only the nonlinear effects up to third order. Let the polarization be split into two 
terms such that 

P{x, t) = Pl{x, t) + Pnl{x, t) (2.8) 

where the linear part and the nonlinear part of the polarization are given as 

PL{x,t) = eo(^j -T)u{x,T)dT^ (2.9) 

Pnl{,xU) = eQ{x^'^^u{xUf) . (2.10) 


is the third order susceptibility and assumed to be a constant. For the linear part, 
since we assume that the medium response is local, retardation in the medium’s re¬ 
sponse to the applied electric field must be taken into account. Note that the relation 


(2.9) is valid in the electric-dipole approximation, which can simply be described as 


the first order approximation. One might expect only the second order nonlinearity 
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to be comparable to the linear response, however, due to inversion symmetry this 
contribution vanishes in optical hbers. We must also note that it is experimentally 
known that the nonlinear effects are weak in silica fibers [18]. For the reason that 


nonlinear effects (third order susceptibility) are small in optical hbers, the nonlinear 


polarization Pnl in (2.8) will be treated as a small perturbation to the total polariza¬ 


tion in the derivation of the wave equation. This leads to a major simplihcation, and 
it is a good treatment for our present disscussion. Therefore, we manipulate equation 


(2.7) with P]\fL = 0 hrst and the include nonlinear term later. Equation (2.7) becomes 


linear when P/vl = 0. If a Fourier transform is applied to the linear equation 


dluix^w) - 1 - —vd^uix^w) = ^o{-w‘^)Pl{x,w) , 


(2,11) 
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then a simple calculation of the Fourier transform of the linear polarization (2.9) gives 


r’+oo 


P]:^{x^w) = / Pi^{x^t)e dt 


/ + 00 ^+oo 

/ “ t)u{x, dtdr 

-oo J —oo 

/ + 00 p-\-oo 

u{x,T)dT / ~dt 

OO J —oo 

^+oo r-\-oo 

= eo / u{x,T)dT / dt 


’ —OO 
/•H-OO 


^H-OO 


= eo / M(x,r)e-™"cir / 

J —OO J —OO 


( 2 . 12 ) 


With this simple form of the linear polarization in hand, we can rewrite equation 


(2.11) as 


dlu{x,w) + ^w'^u{x,w) = ^(— u{x,w) (2-13) 

where the relation c = 1/y/JI^ is used to obtain the coefficient on the right hand side. 
Now the question is how to deal with the linear susceptibility in this equation. 
In general, the linear and nonlinear susceptibilities depend on the frequencies of the 
applied helds. Since there is an assumption of an instantaneous response in the 
nonlinear part of the polarization, we therefore take the nonlinear susceptibility to 
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be a small constant. However, the response of the linear part of the optical held is 
local, and the linear susceptibility is a function of the atomic transition frequency. 
The context in which we are describing propagation of a pulse in a nonlinear cubic 
medium is a semi-classical theory. This indicates that molecules and atoms of the 
medium in which the pulse propagates are explained by quantum mechanics, and light 
itself is governed by the laws of classical electrodynamics. Therefore, the quantum 
mechanical density matrix method applies to the derivation of the linear susceptibility 


in equation (2.13). Because the derivation of the linear susceptibility is tedious and 


requires a good deal of quantum mechanics, we ask the reader to refer to the literature 
EH for the detailed derivation. Instead, we mention briefly how to approximate the 
linear susceptibility. Notice that polarization P is parallel to the electric held E in the 
medium. It is not hard to see as a consequence of this that the linear susceptibility 
can be expressed as a scalar (P = e^xE). Loosely speaking, if the material is modeled 
as a free atom interacting with an electro-magnetic held, the linearized susceptibility 
of the medium in Fourier space can be approximated as 




E/” 


Co 


-2i7naW 


(2.14) 
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where fna is the oscillator strength, cq is some constant, Wna is the resonant frequency 
of the medium and 7„a is a small damping coefficient added to keep the susceptibility 


hnite at the resonant frequency. Although equation (2.14) is a simplified version of 


the linear susceptibility, the exact form of is not important for our results. What 
is important is that we make an approximation of such that can be expressed 
as a polynomial in A. We study the propagation of light in the infrared range with 
wavelengths of 1600-3000 nm, where nm stands for nanometer and the conversion 
relation is Inm = 10“®m. In this regime for silica hbers, a further y^^^ approximation 
can be done such that the linear susceptibility can be approximated to a form of 


yd)(w) y^^^A) = y[,^^ - y^^^A^ 


(2.15) 


where yg^^ = 1.1104 /im ^ and y^^^ = 0.011063 /im Numerical studies based on 


the experimental data show that equation (2.15) is a very good approximation of the 


equation (2.14) for the pulses with wavelengths ranging from 1600 nm to 3000 nm. 


Using the relation between frequency and wavelength A = 21 ^ 0/00 and substituting 
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equation (2.15) in the equation (2.13) leads to 


1 + 

+ ^2 


w’^it — {27 iYx^2^'^ = 0 


(2.16) 


This is the wave equation in the Fourier domain with the linear component of polar¬ 
ization only. As a hnal step, we apply the inverse Fourier transform to the equation 
hrst and add the small perturbation (nonlinear part of the polarization) to the 
equation later. The wave equation in its hnal form ultimately becomes 



dlu = \d^u + \u + 


(2.17) 


where Ci = c/-y/l -I- = 2.065 x 10®m/s, C 2 = = 1.59/im, c = l/^/^o^o 


is the speed of light and is the nonlinear susceptibility. Equation (2.17) is the 


Maxwell equation in one dimension for pulse propagation in the infrared regime 1600- 


3000 nm. 
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2.3 Rescaling the Nonlinear Wave Equation 


In chapter four, we derive an approximate equation replacing the nonlinear wave 
equation for pulse propagation in nonlinear cubic media in the infrared regime with 


the set of equations in (2.6). The new governing equation arises from a scaled version 


of the Maxwell equation. Therefore, we show in this section how to scale the Maxwell 


equation (2.17) in a physical way. 


Let us make a coordinate transformation of the form {x, t) —)■ (^, r) with 


X = x^, t = tr 


(2.18) 


where 
in the 


X and t are properly chosen constants. The double time 


wave equation (2.17) must be modihed according to the 


and space derivatives 
new coordinates such 


that 


_ 1 

dx ‘2 x ‘2 ^^2 
_ 1 
dU 


(2.19) 

( 2 . 20 ) 


Inserting the derivatives (2.19) in the wave equation (2.17) yields an equation in the 
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new coordinates 




( 2 . 21 ) 


From a mathematical point of view, it does not really matter how to choose the 
constants x and t. However, the physics of the problem imposes conditions on these 
constants, and they cannot be chosen arbitrarily from the physics point of view. 


Noting that the nonlinear wave equation (2.17) describes the short pulses propagating 


in a nonlinear cubic medium, these constants must be determined based on short pulse 
parameters. In addition, there is also another condition we must take into account. 
As we will see in chapter four, the coefficient of the Urr term has to be set to one 
in order to derive a short pulse equation over a slow time scale that is introduced 
through a multiple scale expansion. Therefore, we first make the choice of 


y2 

^ 1 P'22) 

as a condition for the derivation of the short pulse equation from the Maxwell equa¬ 
tion. To incorporate the physics of the short pulse equation into the rescaling, let 
us briefly discuss the physical description of a short pulse. A short pulse with wave¬ 
length A = 1.55 /rm has an angular frequency w = 1.94 x 10^^ rad/s { Xu = c). The 
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relationship between angular frequency and the period of a wave can be found using 
w = 2ttv = 27r/T. The period of a short pulse with an angular frequency 1.94 x 10^^ 
rad/s equals to T = 27i/w = 32.4 femtoseconds ( fs ). This means we are dealing 
with a femtosecond time scale in the short pulse dynamics. Hence, we choose the time 
coefficient t to be comparable to the femtosecond time unit such that t = 1 fs. One 
can, on the other hand, choose another t value in the femtosecond unit differing from 


the one we have made here. The choice oii = 1 fs and the condition of (2.22) impose 


the value of x. A basic calculation gives the value of a; as a; = 2.065 x 10 "m ~ 206 


nm. Before computing all the coeffients of the rescaled wave equation (2.21), one 


has to know the physical value of the nonlinear susceptibility. In silica hber, non¬ 
linear susceptibility is given as 1.28 x 10“^® m?/W for a wave with a wavelength 
A = 1.55 jjim [2S1 El]- Using the values i = 1 fs, x = 206 nm and = 1.28 x 10“^® 
/W, one can calculate the unitless coefficents 


^ 2 

a=% = 0.01678573 

b = = 0.005431808 . 

U 


(2.23) 
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Finally, the nonlinear wave equation (2.21) is re-written in its final form after renam¬ 


ing independent variables ^ & r as a; & t, respectively, as 


Uxx = Utt + au + h{u^)tt 


(2.24) 


where a = 0.01678573 and b = 0.005431808. Notice that the coefficient 5 is a measure 
of nonlinearity in this new scale. This is the form of the Maxwell equation we will use 
in the derivation of the deterministic and the stochastic short pulse equation. Before 
proceeding into the next chapter, it should be noted that we shall always refer to this 
equation as the nonlinear wave equation or Maxwell equation in the rest of the work. 



Chapter 3 


Nonlinear Schrodinger Equation 


We discuss the multiple-scale expansion method and the derivation of the nonlinear 
Schrodinger equation (NLSE) in this chapter. The sections contained herein will also 
highlight NLSE solitons and the limitations of the NLSE. 

3.1 Multi-Scale Expansion 

The method of multiple scales is a widely used technique to solve a wide variety of 
problems in physics, engineering and applied sciences [35l |36]. When regular per¬ 
turbation approaches fail, multiple scale expansion is used to approximate periodic 


24 
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solutions to differential equations. In addition to the NLSE and the SPE, the wide 
spectrum 1371 of the application of multiple scaling includes problems such as weakly 
linear and nonlinear vibrations governed by differential equations, the Earth-Moon- 
spaceship problem in the context of orbital mechanics, the role of different time scales 
in flight mechanics, the propagation of waves on a spherical shell in the context of 
solid mechanics, the investigation of the evolution of multi-phase modes for the Klein- 
Gordon equation, the interaction of random waves in dispersive media, propagation 
of nonlinear waves in a cold plasma and approximation for the Thomas-Fermi model 
in atomic physics. For the reason of such a diverse application of the multiple scale 
expansion technique, we shall discuss the method in further detail so as to provide 
gainful insight into the subject. Let us commence with the classical example to 
illustrate multiple scales. We follow the standard example of a damped harmonic os¬ 
cillator. This can be found in any conventional book regarding pertubative techniques 
such as ’’Introduction to Perturbation Methods,” by M.H. Holmes |38]. Consider 


y + ey + y = li, for t > 0 , 


(3.1) 
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where the initial conditions are given as 

?/(0)=0 and 2/'(0) = 1. (3.2) 


This is the equation of a damped oscillator whose analytical solution is well known. 


The exact solution of (3.1) is given as 


^ ^ 2/a ^ sin(Vl-eV4) (3.3) 

Vl - e /4 

Notice that e is a small parameter. We will get the solution of the above equation by 
both the regular expansion technique and multiple-scale expansion, and then compare 
them with the exact solution. Such a comparison illuminates power and beauty of 
the multiple scaling. Let us assume that the function y(t) can be expanded as 


!/(*) = yo(t) + + £^ 2(0 + ■■■ 


(3.4) 
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If equation (3.4) is inserted into the equation (3.1) and the terms of different orders 


in e are rearranged, we will hnd 


y'oit) + yo{t) + e[yl{t) + + y[{t) + y 2 {t)] + ... = 0 (3.5) 


To satisfy the equality, the terms of each order of e including the zeroth order must 
be set to zero. In the zeroth order, we will arrive at the following equation 


y'ait) + yo(t) = 0 


(3,6) 


with the initial conditions ?/o(0) = 0 and 2/o(0) = 0- This is the equation of very 
well-known case, namely, classical undamped simple harmonic oscillator with a unit 


frequency. Equation (3.6) is analytically solvable and the exact solution is found for 


the given initial conditions to be 


yo(t) = sin((). 


(3.7) 
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In a similar fashion, setting the term of the order e to zero gives 


!/i(«) + = -vT) 


(3.8) 


with the initial conditions 2/i(0) = 0 and 2/i(0) = 0. This represents the equation of 
the damped simple harmonic oscillator, and its form is same as that of the equation 


(3.1). The damped harmonic oscillator can be solved by looking for trial oscillatory 


solutions of the form exp(f) or Acos{t) + Bsm(t) because these functions reproduce 
themselves when differentiated. Using the trial form of a A cos(f) + B sin(t) solution 


along with the initial conditions generates the exact solution of the equation (3.8) as 


Vi(t) = --tsin{t). 


(3.9) 


Noting that e is a small constant, we truncate the expansion (3.4) at the hrst order 


and write the approximate solution for the equation (3.1) as 


v(t) “ yi>(t) + = sin(«) - -€(sm((). 


(3.10) 
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In comparison with the analytical solution (3.3), we observe a significant difference 


between the two. The regular expansion approximate solution (3.10) grows in an 


unbounded manner in time because of the second term. The factor et in the second 


term on the right-hand side of equation (3.10) can take values such that it can make 


the whole term as large as the first term or even larger. For this reason, this term is 
called the secular term. On the other hand, the analytical solution decays exponen¬ 
tially with time. The condition required for having a valid approximate solution via 
the regular expansion method is f -C 1/e. 

Let us now solve the same problem using multiple scales. The exact solution 


(3.3) has an oscillatory part that occurs in a time scale on the order of e*^, i.e., 0(1). 
However, the exponentially decaying part has a dependence of et so that it takes 
place on a time scale of the order of 0(l/e). This means the envelope of the exact 
solution, which is governed by the exponential function, changes much more slowly 
than the oscillatory part of the solution, which is governed by the sine term. One can 
then articulate the need for two different scales in the expansion instead of one. This 


is precisely what we are embarking upon now. Let us introduce two time scales such 
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that 


ti = t (3.11) 

t 2 = et (3.12) 


Notice that since we truncate the regular expansion (3.4) at the second term, we keep 


the multiple scaling at two time scales only. In the general case, the exponential term 
depends on et, e^t, as it can easily be seen by expanding the exponent in the 

series form. Therefore, the general case requires multiple scales in the form of 
Note also that the new variables ti and t 2 will be treated as independent variables. 


As a consequence of introducing new variables (3.11), derivatives in equation (3.1) 


must be changed according to the chain rule so that the hrst derivative now becomes 


dti d dt2 d 

dt dti dt dt2 

A A 

dti ^ dt2 


Ft 


(3.13) 
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and the second derivative becomes 


— TTVT + 2e 


d 


+ e 


dU dtf dtidt2 dtl 


(3.14) 


where dtijdt = 1 and dt 2 /dt = e are used. Once we insert (3.13) and (3.14) into the 


equation (3.1), we obtain 


1 d" 


— 

dti dtidt2 


9 

'd^ d' 

+ 6^ 

dti ^^2 


y{t) = 0 (3.15) 


where the initial conditions are y{t) = 0 and {d/dti + e d/dt 2 )y{t) = 1 for ti = ^2 = 0 


( see the derivative relation (3.13) ). Let’s now replace the expansion (3.4) with the 
form 

yif) = yoitiM) + ei/i(^i,0) + e^ 2 / 2 (ti,t 2 ) + ••• (3.16) 


Inserting the equation (3.16) into the equation (3.15) and rearranging the terms lead 


to an equation at the zeroth order of e 

d‘^yo{tiU2) 

dt\ 


+ yoiUi h) — 0 


(3.17) 
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with initial conditions i/q = 0 and dyo/dti = 1 at = ^2 = 0. Notice that we have a 
partial differential equation now, whereas we previously had an ordinary differential 
equation in the regular expansion method at the zeroth order. Partial differential 
equation formalism will actually prevent the secular term appearing in the solution 


as it will be seen shortly. The equation (3.17) is the equation of a simple harmonic 


oscillator as before. The general solution for this equation can be expressed as 


yoih, t 2 ) = 00 (^ 2 ) sin(ti) + 60 (^ 2 ) cos(ti) 


(3.18) 


where ao(0) = 1 and &o(0) = 0- The coefficients 09 (^ 2 ) and 60 (^ 2 ) are functions of ^2 
and will be determined after we handle the next order. At the next order, i.e., 0(e), 
we have the equation 


dtj 



2/1 (^1,0) 


A 

‘dtidt2 dti 


I/o(OT2) 


(3.19) 


where the initial conditions are yi = 0 and dyi/dti = —dyo/dt 2 at 0 = ^2 = 0. If we 
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substitute equation (3.18) in the equation (3.19), we find 


d^yi 

dt\ 


+ yi = ( 260 (^ 2 ) + &o) sin(ti) - ( 200 (^ 2 ) + 00 (^ 2 )) cos(ti 


(3.20) 


The solution of this equation is given as 


2/1 = 01 ( 2 : 2 ) sin(/:i)+ 61 ( 2 : 2 ) cos(ti) 

- ^( 2 &o(^ 2 ) + ho)ti cos(ti) - ^( 200 (^ 2 ) + ao(^ 2 ))ti sin(ti) (3.21) 


where ai(0) = 0 and &i(0) = 0. Remember that we had a secular term in the first 
order solution of the regular expansion. Therefore, there is no reason for there not to 


be one or more secular terms in the solution (3.21) as well. One can observe the direct 


dependence of the solution (3.21) on ti and argue that the terms with ticos{ti) and 


ti sin(2:i) in (3.21) are secular terms in an analogous way to the treatment of the term 


fsin( 2 :) as a secular term in the solution (3.10). Removing these terms will prevent 


the growth of the first order solution (3.21). The secular terms can be disregarded 
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by choosing 


260(^2) + ^0(^2) — 0 


(3.22) 


200(^2) + 00(^2) = 0 

It is not difficult to find the solution of these two equations. Using the initial condi¬ 
tions ao(0) = 1 and 60(0) = 0, we obtain the solutions as 


aoih) = e 

boih) = 0 


(3.23) 


This is the solvability condition of the first order. This means one has to choose 


the coefficients as in (3.23) in order to remove the secular terms appearing in the 


solution (3.21). Notice that the secular terms can be removed without solving the 


equation (3.20) by just setting the coefficients of the trigonometric functions on the 


right-hand side to zero. If we truncate the series (3.16) at the first order and tie 


everthing together, the approximate solution of the problem (3.1) obtained by the 
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multiple scale expansion in the order of e will be 

y{t) Ri yo{ti,t 2 ) + eyi{h,t 2 ) 

= e“"*"/^sin(ti) + e|/i(ti,t2) 

^ e""*"/^sin(ti). (3.24) 


When this solution and the exact solution are compared, one can see that the solution 


(3.24) is a much better solution than the one expressed in (3.10). Note that for 


solutions up to the order of e^, one has to introduce one more time scale, namely, 
ta = The introduction of one more time scale would give an approximate solution 
that is valid over a larger time scale (that is, up to 0(l/e^)), but not necessarily more 
precise over time interval up to 0{l/e) than the approximate solution with two time 


scales. We must also note that two time scales U and t 2 in the solution (3.24) are 
called fast scale and slow scale respectively. 

In general, the time scales depend on the nature of the problem. More complex 
time scaling may be required for some nonlinear problems. Although the exact form 


of different scales may not be clear immediately, multiple scale expansion is a powerful 


method to remove the secular terms in the expansion of the approximate solutions of 
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many nonlinear problems. 


3.2 Derivation of the NLSE 


The nonlinear Schrodinger equation is a model that describes physical nonlinear sys¬ 
tems and may be applied to many nonlinear phenomena such as nonlinear optics, 
nonlinear acoustics, quantum condensates, hydrodynamics and heat pulses in solids 
HI. There are many different forms of the nonlinear Schrodinger equation available 


in the literature. The nonlinear Schrodinger equation has the general form 



(3.25) 


where 2 ; is the propagation direction, i = \/—l, f is the general function that rep¬ 


resents the problem, is the Laplacian operator that can be in one, two or three 


transverse dimensions. Here, we study the NLSE in the context of propagation of 
pulses in a cubic nonlinear medium in one transverse direction. It is a very well-known 
fact that the propagation of optical pulses in optical fibers is governed by the NLSE 
|18j . In the case of a cubic nonlinear medium such as optical fibers, the function / is 
given as f{\u\‘^) = ylup, where 7 is a constant coefficient. 
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One can derive the NLSE by using multiple scales [TTl [13] . If the multiple scale 
expansion of the form 


u{z, t) = eg{z, t) + eg*{z, t) + ... 


g{z,t) = Ao{zuZ2,...,ti)U^^^°-^^°'> 


(3.26) 


Zn = e'^z, to = t, ti = et 


is inserted in the Maxwell equation (2.24), we obtain the nonlinear Schodinger equa¬ 


tion up to a scale transformation in 1 -|- 1 dimension as 



1 

2W ^ 



(3.27) 


where u{z,t) is the magnitude of the applied held (electric held) and is a complex 
quantity, z is the propagation direction, and t is time. The choice of sign depends 
on the cubic (Kerr) coefficient of the nonlinear material. The one dimensional NLSE 
describes the wave propagation in huids and plasmas as well and emerges as a mean 
held equation for many-body bosonic systems in quantum held theory 


The hrst term in the right hand side of the equation (3.27) describes the ehects of 


dispersion. When acting by itself, it does nothing to change the frequency spectrum 
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of the pulse and serves only to broaden or narrow the pulse in time. However, when 
the nonlinear term (the pulse intensity envelope times u itself) acts by itself, it does 
nothing to change the pulse shape in time and serves only to broaden or narrow the 
pulse in the frequency domain. 


3.3 The NLSE Solitons 


The nonlinear Schrodinger equation is a simple partial differential equation with com¬ 
plete integrability |39]. Therefore, NLSE can be solved exactly using the inverse scat¬ 
tering method [ 10 ]. Exact solutions of the NLSE are solitary wave type. This means 
the effects of dispersion and nonlinearity cancel one another and this balance prevents 
the pulse from broadening or blowing up as it propagates [ID] . Solitary wave solutions 


are also called solitons. The choice of the sign in the equation (3.27) gives rise to the 


different type of solitary waves. The soliton solution corresponding to the choice of 
-|- sign is called a bright spatial soliton, whereas the soliton solution corresponding to 
the choice of — sign is called a dark soliton. These types of solutions are sometimes 


called spatial solitons. 
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Let’s consider (3.27) with the plus sign for the nonlinear term; 


du 1 d‘^u . ,2 


(3.28) 


This equation is exactly integrable [39] and the analytical solution in the general form 
na is given as 


u{z,t) = asech[a{t — vz)]exp[ivt + i{a^ — v‘^)z/2] (3.29) 

where a is the soliton amplitude and v is the velocity of the soliton. This is the bright 
soliton solution characterized by the two parameters a and v. The general solution 
reduces to a special solution known as the fundamental solution in the limit v = 0. 
The fundamental solution in this limit would be 

u{z,t) = asech{at)exTp{ia‘^z/2). (3.30) 

This is a particular pulse of an amplitude a whose mean frequency is the central 
frequency. The hyperbolic secant function determines the shape of the envelope 
of the solitary wave, and the exponential function is the phase term. The phase 
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term in the solution (3.30) has no dependence on t, and therefore, the soliton is 


completely nondispersive. This means that its shape does not change with 2 ; either 
in the temporal domain or in the frequency domain [13]. The dispersive term, which 
affects the pulse only in the time domain, and the nonlinear term, which affects the 
pulse only in the frequency domain, cancel each other in a way that stable propagation 
of solitons occur just leaving only a phase shift of the whole pulse behind. The 
fundamental soliton represents the fundamental mode of the optical waveguide created 
by the propagating pulse. If the input pulse has the correct shape in a way that it 


satishes the relation (3.30), all of its energy will be contained in this mode, and the 


pulse will propagate without change in its shape. 


If the sign of the nonlinearity is taken to be minus in (3.27), the nonlinear 
Schrodinger equation becomes 


. du 1 d'^u 


u \ u = 0. 


(3.31) 


This equation is also integrable [ 13 ] and gives rise to the solution in the general form 


u{z,t) = Uo[Bta.iah.{uoB(t — Auoz)) + iA]exp(—itigz), 


(3.32) 
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where the parameters A and B satisfy A"^ + = 1 and Uq is the background am¬ 


plitude. This solution is called a dark soliton. To write the solution (3.32) in terms 


of one parameter instead of two, the relations A = sin cj) and B = cos (j) can be used. 
Here the angle 0 is the half the angle of the total phase shift of 20. When the dark 
soliton is characterized by 0 only, we can write the solution as 


|m|^ = Mo[l — cos^ 0 sech^(Mo cos 4>{t — Uq sin 02 :))]. 


(3.33) 


The term uq sin 0 represents the velocity of the soliton in the z direction. The special 
dark soliton solution can be obtained by setting 0 = 0 such that 

u{z,t) = Mo tanh(Mot) exp(—W q^). (3.34) 


This special dark soliton does not move against the background. Therefore, it is a 
stationary dark soliton and given the name, black soliton. The reason why it is called 


a black soliton is because the soliton (3.34) has a vr phase shift at a; = 0 and the 


intensity drops to zero. In the general case (3.32), the intensity does not drop to zero 


at the center and the soliton is refered to as a gray soliton in such situtations. The 
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phase of a dark soliton changes across its width unlike the phase of a special bright 
soliton which remains constant during propagation. 

Finally, we should note that in many situations the physical constants such as the 
speed of light and pulse width are incorporated in the dehnition of special soliton 
units, which measure distance, time and power only na. The special soliton units 
make nonlinear equations look simple. As a consequence of this simple manifestation, 
adaptation of numerical schemes and understanding solitonic features become easier. 
There are also other types of NLSE solitons such as dispersion managed solitons and 
temporal solitons. The reader can End many works related to the other types of 
NLSE solitons and their interaction [anE]. 

3.4 Limitations of the NLSE 

In this section, we discuss conditions under which the nonlinear Schrodinger equation 
fails to be an approximation to Maxwell’s equations. The derivation of the NLSE 
is based on the multiple scaling technique as shown before. In the derivation of 
the NLSE, we make the assumption that the nonlinearity of the polarization has a 
small contribution to the total polarization and is treated as a perturbation in the 
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expansion of the total polarization. There are cases in which the conpling of the 
freqnencies generated by the medinm may increase the effects of the nonlinearity on 
the propagation of the pnlse, and phase mismatch may occnr as a resnlt BD]. In snch 
cases, the nonlinear Schrodinger eqnation fails. The details of the breakdown of this 
assnmption is beyond the scope of this work. 

The NLSE is a scalar eqnation. To derive the scalar NLSE, we also make the 
assnmption that the electric held maintains its polarization along the hber length. 
This is a fairly good approximation and is valid in many practical sitnations m- 

The third assnmption that is made in the derivation of the NLSE is the slowly 
varying envelope approximation (SVEA) in the propagation direction. This approx¬ 
imation is sometimes called paraxial approximation. The paraxial approximation is 
valid if the width of the pulse is much longer than the light (pulse) wavelength |TT] . 
The slowly varying approximation separates the rapidly varying part of the electric 
held from the slowly varying envelope. Such a seperation is not possible in deriving 
the NLSE if the envelope is assumed not to vary with propagation direction z on a. 
scale much longer than the wavelength. As the pulse width decreases, this approx¬ 
imation begins to breakdown, and the NLSE is not the governing equation of the 
ultra-short pulses anymore [Ml [12]. In the present discussion, the third assumption is 
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the cause of the breakdown of the nonlinear Schrodinger equation in describing the of 
ultra-short pulses. In the next chapter, we will derive a governing equation to replace 
the NLSE for ultra-short pulse propagation in a nonlinear cubic media. 



Chapter 4 


Short Pulse Equation 


This chapter is one of the core chapters of this thesis. We will hrst show the derivation 
of the short pulse equation (SPE). The analytical solution of the short pulse equation 
will be discussed next. We will complete this chapter with the discussion of numerical 
analyses and higher order terms of the short pulse equation. 
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4.1 Derivation of the SPE 

In the derivation of the short pulse equation, we begin with a multiple scale expansion 
ansatz of the form [12] 


u{x,t) = eAo{(j),xi,X 2 ,...) + e‘^Ai{(j),xi,X 2 ,...) + ... 


(4.1) 


with 


t — X 


x„ = e X. 


(4.2) 


The idea is to observe the effects of dispersion and nonlinearity on the pulse over a 
different time scale and see if there exists an equation that is easier in the mathemat¬ 
ical sense and that can model the pulse propagation in the new scales. In this section 
we only consider the expansion up to the order of e, which is a small constant expan¬ 
sion parameter. This means we only introduce two space scales Xq = x and Xi = ex. 


Note that we treat these new variables independently. The equation (2.24) includes 


double time and space derivatives of the u function with respect to the original scales 
X and t. These derivatives must be changed with respect to the new variables in 


the application of the multiple scale expansion to the Maxwell equation (2.24). The 
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transformation of the derivatives follow the chain rnle once again, and the space and 
time derivatives are written respectively as 


d d(p d ^ dxi d 1 d ^ d 

dx dx d(j) dx dxi e d(p dxi 

d d(j) d dxi d Id 

dt dt d(f) dt dxi e dcf) 


(4.3) 


where dcp/dx = —1/e, dxi/dx = e, d(j)/dt = 1/e and dxi/dt = 0, and the second 
derivatives becomes 


d^ 1 d^ 


dx'^ 

dU 


- 2 - 


d d 


e^ dcj)'^ d(j) dx 
1 d^ 


+ C 


dx\ 


(4.4) 


d(jT 


If the relations (4.1) and (4.4) are snbstitnted into (2.24), we obtain 


{eAo + e^Ai + ...) - 2 — — (e^o + + ...) 


e^ d(/)^ 


+e^ 


dxf 


d(p dx\ 

1 ( 9 ^ 

(ev4o + e^Ai + ...) = (e^o + + •••) 


(4.5) 


b d^ 


+a {eAo + e^Ai + ...) + (e^o + + ...) . 


Once we rearrange all the terms, we see that the terms np to 0(l/e) and 0(1) cancel 
ont. As mentioned before, we only keep the terms np to 0(e). Collecting the terms 
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of 0(e) on each side results in 


d , 

e I „ ~ — ^0 


d(/)^ d(p dx 


1 




( 4 . 6 ) 


Finally, if we cancel out like terms, we obtain the following equation 




( 4 . 7 ) 


where the coefficients a and b remain the same with a ~ 0.0168 and b ^ 0.00543. 
This equation is called the short pulse equation (SPE) and derived by T. Schafer and 
C. Wayne [12]. The SPE rests at the heart of the ultra-short pulse dynamics. It is 
the governing equation for ultra-short pulses in silica hbers in the infrared regime. 


4.2 Transformation of the SPE 


We will now show how the original form of SPE (4.7) can be transformed to the form 
of 

UxT = U + -{U^)xx 
6 


(4.8) 
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We start with a general discussion of transformations of differential equations 
and bring the transformation of SPE into this discussion along the way. The exact 
solution of the SPE has been derived simply for this form IZD]. The exact solitary 


wave solution cannot be accepted as the exact solution of the original SPE (4.7) 


without making a proper transformation between the two forms. The transformation 
of SPE is of great importance especially when we check numerically whether the exact 
solution satishes Maxwell equation from which we derive the SPE through multi-scale 
expansion. 

The system of equations 


Mil = ((ag(n) + ± g'(u) 


(4,9) 


»"(«) + w(«) = s 


(4,10) 


are given by M.L. Rabelo to describe pseudo-spherical surfaces. Here a, /3, /x 
and 6 are arbitrary constants, and g(u) is a solution of the linear ordinary differential 


equation of (4.10) and g'(u) stands for the derivative of g with respect to u. The 


partial differential equation of the form (4.9) can be transformed to the form (4.8) 
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for U{X, T) through the following transformations 


u = CiU{X, T) + C 2 , X = c^x + c^t, T = c^t 


( 4 . 11 ) 


where Ci,C 2 ,...,C 5 are properly chosen constants with C 1 C 3 C 5 7 ^ 0 and a 7 ^ 0 . 


First, we will determine whether we can write the original form of the SPE (4.7) 


in the form given by (4.9). If /i = 0 in equation (4.10), then we have 


g{u) = + 'yu + 6 


( 4 . 12 ) 


with arbitrary constants 7 and 6. Once we substitute g{u) and its hrst and second 


derivatives with respect to u in equation (4.9) for which we also choose the positive 


sign in the last term, we obtain 


Uxt = ^ (uxi^aOu^ + ayu + a6 + fl)] + 0 m + 7 . 


( 4 . 13 ) 


Let’s further differentiate each term in the parenthesis. After collecting like terms. 
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we find 


^xt ^xx 


-aUu + aju + 




( 4 . 14 ) 


Choosing arbitrary constants 7 = 0 and ad + (I 


0 rednces eqnation (4.14) to 


Uxt = Ou + aOuiuxY + -aOu^u^ 


( 4 . 15 ) 


The short pulse equation and equation (4.15) are actually equivalent. A simple ma¬ 


nipulation of the SPE is required to show this equivalency. In carrying this out, the 


last term of (4.7) can be expanded in the following manner 


~ ^?i[3(Ao)^(Ao)0] — 6{Ao){Ao)tj)^ + 3(Ao)^(Ao)</,</,, (4.16) 


and the short pulse equation (4.7) can be re-written as 


(Ao)</,xi — Aq{A^)'^ + 2 


(4.17) 


By comparing equations (4.17) and (4.15), one can easily see that these equations are 
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equivalent if we have 


u = Aq, X = <P, t = xi, 6 


a 3b 6b 

2 ’ 


( 4 . 18 ) 


This concludes the transformation of equation (4.7) to the form in (4.9). We will now 


cast equation (4.7) into the form given by (4.8). To accomplish this, we choose the 


coefficents in (4.11) to be 


1 



C2 = 0, C3 


1, C4 = 0, C 5 = 6 


(4.19) 


and rewrite equation (4.11) as 


1 



U{X,T), 


X = x, 


T = et. 


(4.20) 


Combining equations (4.18) and (4.20) leads to the transformation of the kind 


Ao = ^U{X,T) = .[^U{X,T), 

Ja V 6b 


= X, x, = It = --T. (4.21) 

u a 


This is the transformation rule needed to shift from one form of the short pulse 
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equation to the other form. It is rather straightforward to check whether or not we 


have obtained the right transformation rule. Once we substitute equation (4.21) into 


equation (4.17) and carry out the derivatives, we obtain the simplihed equation 


= U+U6U(Uxy + 3U^Uxx) . 


(4.22) 


The second term on the right-hand side is the same expansion we observed in equation 


(4.16). As a hnal step, we write equation (4.22) in its simplest form 


UxT = U+-U^ 


XX 


(4.23) 


To summarize, the original form of the SPE may be transformed to a new form 
through a transformation of the kind such that 


—2d^dxj^AQ — qAq + 6c1|(Aq) 


UxT = U+-U^xx 
6 


(4.24) 


Ao{(j),xi) = 


^U{X,T), 0 = X, 


Xi = 


--T. 

a 
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4.3 Analytical Solution of the SPE 


The short pulse equation (4.8) is an integrable equation [19]. Therefore, its analytical 
solution can be derived by several different methods. We do not show in this section 
how the derivation of the exact solution can be performed. However, we give a brief 
discussion and interpretation of the analytical result. The proof of integrability shows 
that the SPE possesses a Lax pair of the Wadati-Konno-Ichikawa type [121 . Such a 
property allows for the method of inverse scattering to be used in solving nonlinear 
differential equations including the SPE. However, the exact solution of the SPE has 
not yet been found using the inverse scattering theorem. The solution has been found 


using a transformation between the short pulse equation (4.8) and the sine-Gordon 


equation [20|. The sine-Gordon equation is a very well known and very well studied 
equation whose analytical solution is a soliton [lamn. The sine-Gordon equation 


=sm(z) 


(4.25) 




CHAPTER 4. SHORT PULSE EQUATION 


55 


is equivalent to the SPE through a chain of transformations [20] of the form 


v(x,t) = {ul + T) 

X = w{y, t), v{x, t) = Wy{y, t) (4-26) 

z{y,t) = arccos(taj^), 


where the subscripts denote the derivatives with respect to x and y. The trans¬ 


formation (4.26) is used to study the short pulse equation. The exact single-valued, 


non-singular solitary wave solution of the SPE in (4.8) has been derived from a bound 


state kink and a anti-kink solution of the sine-Gordon equation (SGE) by means of a 


transformation of the kind in (4.26), and is given as 


U 

X 


Amn 


m sin -0 sinh cf) + n cos 'ip cosh (j) 


Y -|- 2mn 


m? sin^ p) + rp cosh^ cp 
m sin 2ijj — n sinh 2(f) 


rrP sin^ p) T rP cosh^ (f> 


(4.27) 


with 


0 = m (F-|-T), 'ip = n(Y — T), n = \/(l — m^), 0<m<l. 


(4.28) 
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The exact solution (4.27) represents a nonsingular pulse if and only if m < nicr, where 


TT 

nicr = sin — 0.383 . 

8 


(4.29) 


m is a parameter that determines how short the pulse is. When m reaches its crit¬ 
ical value, the pulse becomes as short as approximately three cycles of its central 


frequency. At the critical value of m, the SPE solution (4.27) becomes singular even 


though it remains single-valued. 

If m is small, for instance m = 0.05, then the solitary wave solution of the SPE 
can be approximated to a simpler form, which can be easily done. As m —)• 0, 
n = VT — m? 1, and when the Taylor series expansion is applied to the hyperbolic 


function sinh0 in (4.27), all sine and sinh terms go to zero. As a result, the solution 


(4.27) can be simplihed to the form 


U ~ 4m 


A F 


cos(F - T) 
cosh(m(F -I- T)) 


= 4m cos(y — T)sech(m(F -|- T)) 


(4.30) 


Notice that the approximate SPE solution (4.30) is similar to the bright soliton of 


the NLSE (3.29) or (3.30). In either case, whether it is the true solution (4.27) 
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or an approximate solution (4.30), the hyperbolic function nonetheless determines 


the envelope of the wave packet while the trigonometric function is responsible for 
oscillations. It should be noted that we will use the single solitary wave solution 


(4.27) in the numerical analysis of short pulse dynamics. 


Multi-soliton solutions of the SPE have also been constructed IZH through a sys¬ 
tematic procedure. This procedure starts with a hodographic transformation rule 
between the SPE and the sine-Gordon equation. The hodograph transformation is a 
method used to transform nonlinear partial differential equations into linear partial 
differential equations. Using the breather solutions (solitons) of the sine-Gordon equa¬ 
tion, the solutions of the SPE can be written as a system of linear partial diffential 
equations governing the inverse mapping to the original coordinates. The analytical 
integration can be done for this linear system of partial differential equations, which 
yields the analytical multi-loop and multi-breather solutions. The one-breather solu¬ 
tion (solitary wave solution) is given by 


u(y,t) = 2i 
x(y,t) = y-2 + d 


(4.31) 
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with 


/ = 1 + + ( “ ] 


f =r 


0 = a l y + 


^i = 0 + iX 


t]+X 


(4.32) 


+ 6 ^ 


X = b(y- 


t ] + y, 


+ 6 ^ 

where a and b are positive constants, y and A are real constants, and d is the inte¬ 
gration constant. Note that the choices of different values of d make a difference only 


in the location of the pulse. Equations (4.31) and (4.32) represent the compact form 


of the one-soliton solution of the short pulse equation (4.8) obtained by a different 


technique than the one used to arrive at the solution of the form (4.27). In order to 


show that the solutions (4.27) and (4.31) are equivalent, we will explicitly write the 


solution (4.31) in the parametric form as [2T 


u{y,t) = 


Aab b sin(x) cosh (0 -|- In — a cos(x) sinh (6^ -|- In 


+ &2 


6^ cosh^ (0 -I- In -I- cosh^(x) 


x{y,t) = y- 


2ab a sin(2x) -|- b sinh (26 -|- 2 In 4a 


(4.33) 


-|- 6^ 62 cosh^ (0 -|- In cosh^(x) a‘^ + 62 


+ d. 


This solution is multi-valued without any conditions imposed. Since we are searching 
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for a single-valued nonsingular solution, the condition giving rise to a one soliton 
single-valued solution is 


/— (2 
— \/2 -|- 1 < — 


cos(x) 


b cosh (6* -|- In 


< 


V2-I. 


(4.34) 


Recall that we have imposed a and b positive constants beforehand. However, the 


condition (4.34) imposes an additional constraint on a and b such that 0 < a/b < 


y/2 — 1. Notice that the phase 6 governs the envelope of the soliton whereas the phase 


X is responsible for the oscillations. The comparison of the solutions (4.27) and (4.33) 


shows that the two are equivalent. The advantage of using the solution (4.33) is that 


we can vary the speed of pulse since the speed is dehned to be c = l/(a^ -|- 6^). This 
flexibility is very useful when studying the solitonic features numerically. 

In a similar fashion, multi-soliton solutions of the SPE can be constructed 


We present the multi-soliton solution here because we intend to discuss the collisons 


of the SPE solitons in chapter Eve. The parametric multi-soliton solution can be 
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expressed in the compact form as 


u(y,t) = 2i 

x(y,t) =y-2 (^ln/7)^ + d 


(4.35) 


with 


N 


f = explj; Tj 2 0 ^ ^ TjTk'yjk] 


fi=0,l j=l 

^<j<N 

N 


f' = {xj - 

1 4“ ^ ^ PjPkXjk] 

1^=0,1 j=l 

^<j<N 

0 = pjy + + Oo, 

Pj 

(j = l,2,...,M) 

= {3,k = 

Pj + Pk 

1,2, ^ k) 

P2j = ttj + ibj, ttj > 0, bj > 0, 

(j = l,2,...,M) 

■ 62 ^- 1,0 = 'C 2 yo = 7 ' 

(j = l,2,...,M) 

9j = aj{y + Cjt) + Xj, 

(j = l,2,...,M) 

Xj ^j{y 4” Pji 

(j = l,2,...,M) 

1 

(j = l,2,...,M), 


(4.36) 


where pj and Qq are arbitrary parameters satisfying the conditions pj 7 for j ^ k, 
i is the imaginary number such that i = V—1, N is an arbitrary positive integer, and 
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M = N/2 is the number determining the multi-soliton solutions (one soliton, two 
solitons, etc). i i^eans summation over all possible combinations of /ii = 0,1, 

/i 2 = 0, l,...,/iAr = 0,1. Notice that when M = 1 (or N = 2), multi-soliton solutions 


given by equations (4.35) and (4.36) reduce to the one-soliton solution in (4.31) and 


(4.32). To derive the two-soliton solution, one must choose iV = 4 and M = 2. The 


condition for a single-valued multi-breather solution is 


M 

0 < ^ ^ < ^2 - 1 (4.37) 

j=i T 

It is worthwhile to note that single and multi-loop solutions of the SPE have been 
derived as well [201 El]- Although these loop solutions show solitonic features, loop 
solutions are not single-valued. We are only interested in the nonsingular, single¬ 
valued solitary wave solutions in the context of, and in application to, nonlinear 
optics. 


4.4 Numerical Analysis of the SPE 

The short pulse equation is an exactly solvable equation whose analytical solution is 


a soliton as mentioned in the previous section. We have numerically shown that these 
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solitary waves persist in the short pulse equation if they are used as initial conditions. 


Figure 4.1 shows that if an ultra-short solitary wave (4.27) is chosen as the initial 



X 


Figure 4.1: The numerical solution of the SPE at t=50 and t=100 units of propagation 
distance. 


condition, this initial wave propagates stably in the SPE (4.8). In this experiment, 


we pick the soliton parameter m = 0.3 so that we have the very short solitary wave. 


For the other values of m such as m = 0.2 and m = 0.05, which accordingly generate 
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solitary waves with different sizes, we observe stable propagation as well. We see the 
initial solitary wave at t = 0 in £gure|4T Note that x is the temporal variable in our 
case. The propagation of the initial pulse is shown at t = 50 and t = 100. 

Let us also compare these results with the results obtained from the exact solution 
to see whether there is any deviation from the exact shape and size of the soliton as 
it moves along the line. By doing so, we compute the exact solutions at t = 50 and 


t = 100 using the analytical solution (4.27), and compare them with the corresponding 


numerical results. 


Figures 4.2 and 4.3 exhibit the comparisons of these two solitons at t = 50 and 


t = 100 respectively. Since there are no observable differences between the exact and 
the numerical results at both t = 50 and t = 100, we also plot the maximum error for 
each case. The error in both cases is almost zero, and therefore we can just relate this 
negligible differences to the experimental error. To reiterate, an initial solitary wave 
obtained from the exact result for any values of the soliton parameter m propagates 
stably in the short pulse equation, and this serves as numerical proof for the analytical 
solution of the SPE. 

An interesting question to pose is to ask what happens if the exact solitary wave 


solution (4.27) is used as an initial condition in the original equation (2.24). It 
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Figure 4.2: Comparison of the exact solitary wave solution to the numerical result at 
t = 50 distance units (upper graph). Maximum error between the exact solution and 
the numerical solution at the same distance (bottom graph) 
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Figure 4.3: Comparison of the exact solitary wave solution to the numerical result 
at t = 100 distance units (upper graph). Maximum error between the exact solution 
and the numerical solution at the same distance (bottom graph) 
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has already been shown |15] that the exact solitary wave of the short pulse equation 


persists in the linear wave equation (6 = 0 in (2.24)). We have also shown numerically 


that these solitons propagate stably in the nonlinear wave equation |46]. Figure 4.4 



Figure 4.4: Evolution of the SPE soliton in the Maxwell equation 


displays the exact initial SPE soliton propagation at a: = 25 and its comparison to the 


analytical result at a; = 25. Note that the evolution variable is switched to t in the 
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nonlinear equation. The solid line shows the analytical solution at x = 25, whereas 
the dashed line shows the numerical solution at the same distance. As one can see 


from hgure 4.4, the initial SPE soliton persists in the nonlinear Maxwell equation. 


It may noteworthy to say a few more words about the details of the numerical ex¬ 
periment. We modify the initial SPE soliton according to the multiple scale expansion 


(4.1) and (4.2) such that the magnitude of the initial pulse becomes u = €u(x/€,0) 
at t = 0. Accordingly, the analytical result has to be modihed as well. Therefore, we 
have the analytical solution at t = 25 distance as Uana = eu ([ (x — 25)/e], —e25). We 
choose the soliton parameter m = 0.3 and the expansion parameter e = 0.2. Note 
also that SPE is the leading order 0(e) approximation and choosing the propagation 
distance 0(l/e^) = 25 units in the numerical experiment is more than enough to 
observe any abnormalities in the propagation. 

Before closing this section, let us also discuss the error accumulation when the 
initial SPE soliton propagates in the nonlinear wave equation. Since the SPE solitons 
are not the exact solutions of the nonlinear wave equation, we expect, in general, a 
growing difference between the numerical and the analytical results with the prop¬ 


agation distance. For instance, the maximum error at t = 25 is 0.0186 unit. A 
more elucidating picture for this discussion is to plot the L^, L? and L°° norms of 
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/ = uNumerical “ ^Exact veisus propagation distance, and the norms of / are defined 
respectively as 


ll/lli- = (/ l/WI*) 


LOO = max(|/(t)|). 


In figure 4.5, the dotted line shows norm, the dashed line shows norm and 


the solid line shows the L°° norm. The solution of the leading order equation of 
the multiple scale expansion fails to approximate the solution of the nonlinear wave 


equation at the larger propagation distances as can be seen by hgure Hug. 


4.5 The Higher Order SPE 


If a multiple scale expansion of the form (4.1) with the scale transformation (4.2) is 


applied to the Maxwell equation (2.24), the evolution of ultra-short pulses in optical 


fibers is expressed by the function y4o(0, Xi) over the scales (0, Xi) instead of the 
function u{x,t) over the scales (t, x). If one wants to improve the accuracy of the 
expansion, one can take into account the dependence of the Ao function on X 3 as 


introduced by the multiple scale expansion (4.1) and (4.2). In the derivation of 









L Norm 
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Figure 4.5: The growth of the deviations for the evolution of the SPE soliton in the 
nonlinear wave equation as dehned by L^, and norms. 
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the short pulse equation, we only consider the expansion up to the order of e, and 
the evolution equation is expressed over two variables, 0 and xi. In this section, 
we consider the higher order expansion terms so that the Aq function has another 
dependence, i.e., Aq[(1),xi) —>■ Ao(0,Xi,Xa). To insure that, we introduce more space 
variables such that X 2 = e^x and xs = e^x according to the scale transformation 


(4.2). It is imperative to note that we follow the same procedure here as we followed 


in section 4.1. Note also that the presence of Xq in the expansion leaves the result 
unchanged. Therefore, we only keep the new space variables xi, X 2 and xs here. If 
the procedure in section 4.1 is repeated, the terms of 0(l/e) and 0(e) canceled out. 


and the terms of 0(e) generate the short pulse equation (4.7). By choosing Ai = 0 
and dAQ/dx 2 = 0 (i.e, Aq is independent of X 2 ), the terms of O(e^) will canceled out 
as well. Finally, the terms of O(e^) become 


a^A„ a^A„ 

d‘^xi d(j)dx3 


(4.39) 


This equation underlines the dependence of Aq on X 3 . Now it remains to be seen how 
exactly one incorporates the X 3 dependence in the short pulse equation. To write a 
single evolution equation for Aq using the terms of 0(e) and O(e^), we will combine 
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the two equations given by (4.7) and (4.39). By doing so, let us first integrate the 


SPE (4.7) with respect to 0, 


(^0)11 



(4.40) 


Notice that we have the double xi derivative in equation (4.39). We now take the 


derivative of (4.40) with respect to xi 



^2 r<i) rcj} ab 3 

— — / / Aod 0 + ^(^ 0 ) + 

^ J —00 J —00 ^ 


(4.41) 
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If we substitue the expression for {Aq)xixi 
it with respect to (j), we obtain 


in (4.41) into equation (4.39) and integrate 


2(^o)x3 







Sab 

~T 


(Aof 



AQd(j) + 




(4.42) 


Introducing a new variable such that m: = xi and e^x = X 3 , we hnd 


— (^o)xi + e^(^o)a;3- 


(4.43) 


Finally, we obtain (Aq)®! and (^ 40 ) 3,3 from the relations (4.40) and (4.42) respectively 


and substitute them into equation (4.43) so that 


(^ 0 ) 





X 0 X 3 

4 



SXoXs .2 
4 ° 


A^dcf) + 



(4.44) 


This represents the higher order short pulse equation. Notice that Ao(^)0) is the 
magnitude of the electric held following the introduction of the new variable x. For 
higher orders, we expect to improve our numerical results. The numerical validation 
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of the higher order SPE remains an open problem, and one in which we hope to tackle 
in the future. 



Chapter 5 


Solitons 


Solitons manifest themselves in many branches of modern science snch as nonlinear 
optics, plasma physics, hydrodynamics and biology They have been nsed 

extensively in optical commnnication [181 US US] since the discovery of the nonlinear 
Schrodinger eqnation. Althongh we have tonched on various soliton solutions in 
the preceding chapters, we provide a historical synopsis of solitons in the present 
chapter, and also delve into a brief discussion of the most famous solitons such as 
KdV solitons. We aim to show and interpret our numerical results and in-so-doing 
validate the solitonic properties of the SPE solitons. 
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5.1 What is a Soliton? 

A soliton is a wave packet or a pulse that maintains its shape as it propagates with 
a constant speed in a medium because of a delicate balance between dispersive ef¬ 
fects and nonlinearity |13]. Although the term soliton was introduced in the 1960s, 
physical solitary waves were hrst observed in water waves by J.S.Russell in 1834. The 
observed nondispersive water waves were just the analog of the latter optical solitons. 
The description of Scott Russell’s solitary water waves was published in his paper in 
1844 [13], which includes the following quote: 

/ was observing the motion of a boat which was rapidly drawn along a narrow 
channel by a pair of horses, when the boat suddenly stopped-not so the mass of water 
in the channel which it had put in motion; it accumulated round the prow of the ves¬ 
sel in a state of violent agitation, then suddenly leaving it behind, rolled forward with 
great velocity, assuming the form of a large solitary elevation, a rounded, smooth and 
well-defined heap of water, which continued its course along the channel apparently 
without change of form or diminution of speed. I followed it on horseback, and over¬ 
took it still rolling on at a rate of some eight or nine miles an hour, preserving its 
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original figure some thirty feet long and a foot to a foot and a half in height. Rs height 
gradually diminished, and after a chase of one or two miles I lost it in the windings of 
the channel. Such, in the month of August 1834, my first chance interview with 
that singular and beautiful phenomenon which I have called the Wave of Translation. 

The name soliton was given later to such wave translation observed in 1834. In 
response to the observation of the wave of translation by John Scott Russell, Boussi- 
nesq’s equation and KdV equation were derived by Joseph Boussinesq, and Diederik 
Korteweg and Gustav de Vries in 1872 and 1895 respectively to describe the wave 
of translation mathematically Ha. Both equations can be solved exactly and the 
simplest solutions turn out to be solitary waves. The third-order KdV serves as a 
model for waves on shallow water surfaces, and can be given in standard form as 

ut + 6uu,, + {u)„^^ = 0, (5.1) 


where z is the propagation direction and t is the time variable. The exact solution of 
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the KdV equation is a soliton and can be written as 

u(z,t) = asech^(\/a/2(z — 2at)), (5.2) 

where a is a constant number representing the amplitude of the initial soliton. 
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Figure 5.1: Comparison of the SPE, NLSE and KdV solitons 
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It may be enlightening to compare all the solitons we have thus far mentioned, 


i.e., the SPE, NLSE and KdV solitons. In hgure |5.1[ we show the SPE, NLSE 
and KdV solitons at the propagation distance t = 50 (in accordance with the SPE 
transformation, t is expressed in units of distance). We directly observe the similarity 


among the KdV soliton (5.2), the NLSE fundamental soliton (3.29) and the SPE 


approximate solution (4.30) in that they all exhibit a similar envelope shape albeit 


the KdV soliton has no oscillatory contribution. It is worthy to note that upon 
varying the parameter m up to a critical value {rricr ~ 0.383), we can alter the width 
of the envelope of the SPE solitons. In reference to the SPE soliton shown in hgure 


5.1, we have chosen the parameter m to be 0.05. With regards to the NLSE soliton 


in (3.29), hgure 5.1 depicts the soliton with a unit velocity and an amplitude of 0.2 


unit. 


The term soliton was introduced by Zabusky and Kruskal in 1965 while working 
on the interaction among KdV solitary waves numerically. These solitary waves were 
named as solitons because of the particle-like behavior when they collide. The exact 
solution of the KdV equation was found by Gardner in 1967 through a method called 
the inverse scattering transform. Other notable nonlinear equations possessing soliton 
solutions, such as the NLSE and the sine-Gordon (sG) equation, may also be solved 
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analytically by the inverse scattering method. 

KdV solitons are different than the optical solitons. They describe the solitary 
wave of a wave. On the other hand, an optical soliton describes the solitary wave of 
an envelope in a nonlinear cnbic medinm. Optical solitons are electromagnetic waves 
that are self-localized or self-trapped. This means they move at a constant speed 
in a medinm with no change in their shape becanse of a delicate balance between 
nonlinearity and dispersion. The hrst technological application of optical solitons was 
done in 1973 for pnlse propagation in optical hbers. Since then, NLSE solitons have 
been nsed in many practical sitnations. Using nltra-short solitons in data transfer 
and commnnication may expand the spectrnm of technology [m [E]. Dne to the 
possibility of a wide-range applicability, we will farther investigate nltra-short solitons 
(the SPE solitons) and show their solitonic properties nnmerically. 


5.2 The SPE Solitons 


An analytical solntion of the SPE is given by (4.27) and (4.30). The latter is an 


approximation of the exact solntion whenever the parameter m takes small valnes. 


The parameter m determines the shape of the pnlse. Fignre 5.2 shows the exact 
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Figure 5.2: The parameter m in the exact one-soliton solution determines the shape 
and width of the SPE solitons. Three exact SPE solitons for m = 0.1, m = 0.2, 
m = 0.35 are shown at t = 10 units of propagation distance. 
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solution (4.27) for different values of m at t = 10. As we increase m from 0 to 0.35 


see (4.29)), the pulse width narrows and the amplitude grows. In the case by which 


m takes on smaller values, the width widens and decreases in amplitude. For higher 
values of m, we observe pulses as short as three cycles of its central frequency. When 


we choose a small value for m, one can use the approximate solution (4.30) instead 


of the exact solution (4.27). The numerical work displayed in figure 5.3 demonstrates 


that this is a very good approximation for m = 0.05. However, if m = 0.1 is chosen. 


the approximation starts to fail as can be seen in hgure 5.3 For this reason, we 


suggest that an approximate solution should be used for m values that are equal to 
or less than 0.05 if necessary. 

As it was emphasized before, the solitary wave solution of the short pulse equation 
lends itself to a stable propagation due to a fine balance between dispersion and 
nonlinearity. Dispersion, on one hand, drives solitons to spread out as they propagate. 
On the other hand, nonlinearity gives rise to a centralizing effect on the solitons 
tending to draw them together. The linear broadening of a soliton is canceled out 
by the nonlinearity of the medium whose origin is the intensity dependence of the 
refractive index [T3]. If one of these balancing effects is lost, the result is an unstable 
soliton which cannot exist over an extended period of time. Our numerical simulations 








CHAPTER 5. SOLITONS 


82 




Figure 5.3: Comparison of the exact analytical and approximate solutions of the SPE 
for m = 0.05 (upper graph) and m = 0.1 (bottom graph) at t = 50 distance units. 
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validate such effects for the SPE solitons. 

Let’s now switch off the nonlinearity in the SPE. Without the nonlinear term, the 
SPE takes the form 

U XT ~ U . ( 5 - 3 ) 


Figure 5.4 displays the broadening of the solitary pulse as it propagates along the 



Figure 5.4: The propagation of the SPE solitary wave in the absence of nonlinearity. 


























































CHAPTER 5. SOLITONS 


84 


medium without the nonlinear term. As anticipated, dispersion dissipates broadens 
the pulse linearly in the absence of nonlinearity. 

Leaving on nonlinearity, and now switching off dispersion, the SPE takes the form 


UxT = -^{U^)xx 
0 


( 5 . 4 ) 


Once we let the pulse propagate in the presence of nonlinearity alone, the nonlinear 
term forces the pulse to be more concentrated at the center as it tends to travel. This 
centralizing effect is not balanced by the linear broadening effect due to dispersion 


and the pulse eventually blows up. Figure 5A shows the centralizing effect of the 
nonlinear term at the propagation distance t = 150. The pulse does not move, and 
it becomes more concentrated in the absence of dispersion. If we allow the solitary 
wave to propagate further, it would blow up. 

We have shown numerically that SPE solitons shows a stable propagation due to 
the key balance between dispersion and nonlinearity. This is a unique property of 
solitons. In the next section, we will simulate the particle-like behaviour of the SPE 


solitons as they collide. 
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5.3 Colliding the SPE Solitons 

Solitons are localized excitations propagating in a system with a constant velocity. 
They behave like particles [13]. When there is a large seperation distance between 
two solitons, they essentially do not interact. Once we allow for the two solitons to 
move in opposite directions towards each other, each moves with a constant shape 
and velocity. As the solitons approach one another, their shapes begin to deform. The 
waves merge into one another. In the process, a wave packet is formed and as such, 
it cannot be represented as a linear combination of two solitons. This wave packet, 
however, soon splits into two solitons each with the same shape and velocity as before. 
Thereafter, the solitons move along in their respective directions as if nothing had 
happened. 

We have numerically validated this property of the SPE solitons by colliding the 
exact solitary waves of the SPE using an exponential time differencing (ETD) method- 
based algorithm. The numerical experiment validating the soliton interaction employs 


the two-soliton solution given by equations (4.35) and (4.36). We choose = 4 to 


derive the two-soliton solution of the SPE. The parameters determining the speed of 
each soliton are chosen such that ai = 0.1, bi = 0.5, 02 = 0.16, and 62 = 0.8. Notice 
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that the condition of having a singular two-soliton solution (4.37) is satished. With 


these particular values of a and b parameters, the speed of the hrst soliton is ci = 1.50 
and the speed of the other is C 2 = 3.85 units. We show the initial two-soliton solution 


in hgure 5.6 The larger pulse is the one with the higher speed. As they propagate 


along the line, the faster soliton collides with the slower one ( the smaller of the two) 
and they merge into each other. They soon split apart into two pulses so that the 


one traveling with the higher speed is followed by the other as shown in hgure 5.6 


The next issue to concern ourselves with is to check whether or not the solitons 
remain the same after they collide. We then compare the solitons with the analytical 
result. 


As we see in hgure 5T, we do not observe any discernable diherence between the 
exact solution and the experimental result. We also show the maximum error between 


these two results in hgure 5.7 since it is almost impossible to distinguish one from the 


other. 


It may also be noteworthy to mention that one can validate the interaction picture 


using the nonlinear wave equation (2.24). In that case, one can use the two one-soliton 


solutions of the form (4.27) as an initial condition. If the sign in the propagation 


variable of the one-soliton solution is made minus, the soliton will be traveling in the 
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Figure 5.6: Two SPE solitons are apart before the collision takes place (upper graph) 
and two SPE solitons after the collision occurs (bottom graph). 
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opposite direction. Therefore, these two one solitary wave solutions can be used as the 
initial condition of the numerical schemes like two solitons moving towards each other. 
We have carried out the numerical experiment in colliding two one-soliton solutions 
by employing the Ablowitz-Ladik algorithm (see chapter seven), however, the results 
obtained are not included herein as they only reproduce the the particle-like property 
of the SPE solitons as previously shown. 



Chapter 6 


Stochastic Short Pulse Equation 


We will derive a stochastic version of the short pulse equation in this chapter. The 
impact of stochasticity on ultra-short pulse propagation and the comparison of the 
stochastic SPE with the stochastic nonlinear wave equation will be discussed via the 
numerical experiments. 

6.1 White Noise and Discrete Noise 

The physical system subject to random fluctuations are modelled by stochastic dif¬ 
ferential equations (SDEs). These stochastic equations may be either in the ordinary 
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differential or partial differential form. In general, the random fluctuations in the 
STDs are repsented by white noise [7], which is a random process. There are, of 
course, other types of random processes such as jump processes. 

Let us consider a simple SDE 

^ = a(x,t)+ , (6.1) 


where is a function of time representing the random fluctuations in the system, 
and a(x, t) and b{x, t) are some given functions depending on the variables x and 
t. This equation is sometimes called the Langevin equation. The function ^{t) rep¬ 
resenting the random fluctuations is called the white noise. The derivative of the 
well-known Brownian motion or Wiener process W{t) may be associated to the white 
noise such that 

or dW{t) = . (6.2) 


By employing the relation between the white noise and the Wiener process (6.2), we 


can write the SDE (6.1) in the form 


dx = a{x, t)dt + b{x, t)dW(t). 


( 6 . 3 ) 
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The treatment of the noise in a physical system can also be analyzed from the perspec¬ 
tive of the Fokker-Planck equation, which is an equation of motion for the probability 
distrubution function [8]. Therefore, the physical systems in which there is fluctu¬ 
ating random noise can be studied through stochastic differential equations or the 
corresponding Fokker-Planck equations. An ideal mathematical formulation of the 
white noise is a stochastic process with zero mean such that 


m) = 0 


(6.4) 


where the Dirac delta function is dehned as 


S(t' -1) 



oo 




(6.5) 


White noise does not exist in the physical world. The noise in our system has actually 
a hnite correlation. The delta function correlation of the white noise is just the 
idealization of realistic noise, and it is indeed a good representation of the physical 
noise in a mathematical sense [7]. 
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It may also be remarkable to emphazise that the Fourier transform of the white 
noise is a diffferent noise, but the noise distribution in the Fourier domain appears to 
be random noise as well 1501. 


The origin of the noise in our physical system will be discussed in the next section. 
We will now show how the white noise as a mathematical representation of the physical 
noise can be discretized. This is particularly important when stochastic equations 
are studied numerically. We assume there exists a discrete approximation y of the 


continuous white noise ( given by (6.4). If the time interval is divided into N intervals, 


we have a discrete random number for each interval. The collection of N random 
numbers xi,X2, ■■■, Xn drawn at times tj = At, t2 = 2At, ...,tN = t has a correlation 
such that 


iXiXj) = i = j = 1,2,...,N , 


( 6 . 6 ) 


where is the variance of the discrete noise and 6ij is the Kronecker delta function 


such that 


6ij < 


0 Hi 
1 if i=j 


(6.7) 


The Kronecker delta function (6.7) is just the discrete analog of the Dirac delta 
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function (6.5). The variable W(t) in equation (6.2) is the continuous Gaussian random 
variable and is the integration of the continous white noise 


W{t) = / at )dt 
Jo 


( 6 , 8 ) 


The corresponding Gaussian random variable can be dehned as 


N 




(6.9) 


Since the discrete noise is the approximation of the continuous noise, the statistical 
properties of the continuous process W{t) and discrete process Y must therefore 
match. The continuous process has a zero mean and a variance of t as given by 


equations (6.4) and (6.8). The mean of the sum of all random variables will be 


assumed to be zero. The variance of the discrete process may not be so evident. 


Hence, 


N 


Var{y) = (At)^ ^ = Na‘^{Aty , 

i=i 


( 6 . 10 ) 


where we have used the relations (6.6) and (6.7). This is the variance of the discrete 
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process and must be the same as the variance of the continuous process. Therefore, 


N(T‘^{Atf = 1 . 


( 6 . 11 ) 


Because NNt = t, the discrete variance 




At 


( 6 . 12 ) 


If we choose the variance 1/At in the discrete case, we will match the statistical 
properties of the continuous process with a variance of t. 

We have to emphasize that the strength of the noise is assumed to be one in this 


discussion (see the relation (6.4)). In the cases where the strength of the noise is not 


one, the variance of the discrete noise would be cr^ = i^/At, where v measures the 
strength of the noise. 


6.2 Derivation of the Stochastic SPE 

The stochastic pictures of the nonlinear models are a fundamental issue in nonlinear 
science. Random fluctuations are widely present in nature and the source of the 
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stochasticity may not be spotted easily. Although the deterministic models may 
work well in many situations, the stochastic perturbations may have non-negligible 
effects in some cases. The propagation of ultra-short pulses in a nonlinear medium 
in which stochasticity is taken into account is a more realistic situation. This is the 
main motivation behind the attempt to find a stochastic model for ultra-short pulse 
dynamics. 

Optical soliton propagation in hbers in the presence of a stochastic perturbation 
has been studied in the context of the NLSE model. The nonlinear Schrodinger 
equation with a linear multiplicative stochastic term is a well known model for pulse 
propagation in nonlinear media that exhibit a stochastic nature. The sources of 
randomness in optical hbers vary. Stochasticity may cause the phase of the wave 
huctuate. The small huctuations in the pulse size or intensity can grow with propa¬ 
gation and may eventually lead to a pulse collapse. The inhomogeneities in a hber’s 
core, or the huctuations in the linear refractive index of the core, are the major 
source of medium-related stochastic phase huctuations. The possible sources of the 
phase fuctuations are stimulated Brillouin scattering, stimulated Raman scattering 
and medium inhomogeneities [SB HD- The stochasticity may come from the non¬ 
linearity of the medium as well. The dynamical ehect of the noise added by the 
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stochastic nature of the nonlinearity is nowhere comparable to the noise due to the 


inhomogeneities 


Nevertheless, there may be other sources of randomness play¬ 


ing a role, and they may originate from the other parts of the system such as the 
inherent power fluctuation in lasers used as input pumps. Apart from these, quantum 
phase fuctuations are also sources of phase noise in optical hbers although they are 
practically negligible [15]. In many aplications, Langevin noise (white noise) is used 
to study the fluctuations in the sytem 


The deterministic equation (2.24) leads to a model that describes ultra-short pulse 


dynamics in a deterministic way (equation (4.8)). To derive a stochastic modeling 


equation for our system, one must hrst consider whether there exists any fluctuations 
in the system. As discussed previously, the main source of the noise appears in linear 
polarization of the medium, and fluctuations in nonlinearity are quite small if one 
considers the pulse propagation in the context of NLSE. Recall that the nonlinear 
part of the polarization is treated as a perturbation to the total polarization. The 
fluctuations in the perturbed term will be ignored here for the reason that already 
a small noise in a small nonlinear term does not play a signihcant role in the pulse 
dynamics. We claim the randomness in nonlinearity is much smaller than the noise 
in the dispersion term as in the case of NLSE and is insignihcant. We can now 
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argue that the linear polarization of the material in response to an applied electric 
held is not the same everywhere in a nonlinear medium, but huctuating. Small 


huctuations then appear in the dispersion term in equation (2.24) allowing us to cast 


the deterministic Maxwell equation into a stochastic form. More rigorously, one can 
add a huctuating term to the approximate value of the linear polarization in the 
Fourier domain. Let us introduce the noise in the dispersion term such that the 


huctuating linear susceptibility in equation (2.15) appears to be 


^ - (X2’ + I^'Xrand)X‘^ , 


= yd) _ tyd) 


(6.13) 


where Xrand represents the small noise in the Fourier domain and 12 ' is the strength 
of the noise. Once we substitute the huctuating linear susceptibility into equation 


(2.13) and follow the same procedure we used to obtain the deterministic Maxwell 


equation (2.24), the rescaled stochastic Maxwell equation can be written as 


Uxx = utt + (a + + h{u^)tt , 


(6.14) 


where v is the rescaled strength of the noise. This is the stochastic version of the 
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nonlinear wave equation. The noise in the system is modelled as a white noise whose 
statistical properties are defined as 


('^(^)) = 0 

= S{x-x') 


(6.15) 


Since the nature of noise in many physical systems exhibits a normal distribution with 
zero mean, white noise offers a convenient mathematical implementation as under¬ 


stood by (6.16). Furthermore, the fact that fluctuations in the linear polarization are 


small and the average over these fluctuations is zero, averaging the stochastic wave 


equation (6.14) removes fluctuations from the system leading to the deterministic 


Maxwell equation (2.24) 


We have now reached a point at which we can derive a stochastic short pulse 
equation. To make a distinction between deterministic and stochastic equations, we 
will make a notational change and replace the amplitude of the applied field u{x,t) 
with E(x,t). The stochastic wave equation, in terms of our new notation, is now 
expressed as 

L/xx = Ett + (n + v^{x))E -|- b(E^)u ■ (6.16) 
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One can often obtain a more useful mathematical expression for a given system by 
introducing a fast scale as well as a slow scale. We now proceed to derive a stochastic 
SPE in a similar manner to the way in which we derived the deterministic SPE in 
chapter four. A multi-scale expansion of the form 


E(x,t) = eM„((t>,xo, xi,X 2 ,...) + x„,xi,X 2 , ■■■) + ■■ 


(6.17) 


with new scales 


t — X 


Xr, = e X 


(6.18) 


is used in the derivation of the stochastic SPE. Notice that Mo(0, xq, xi, X 2 ,...) has 


a dependence on xq, whereas the Ao(0, Xi, X 2 ,...) term in expansion (4.1) for the 
deterministic case has no dependence on the scale xq. If xq had been incorporated 
in the Aq function (Aq = Ao(0, Xq, Xi, X 2 ,...)), the same deterministic equation (the 


SPE) (4.7) would have been derived. The reader’s curiosity may demand to know 


what occurs if xq is removed from the expansion (6.17) as in the deterministic case 


see 


4.1). We will soon see that this is not permissible in the stochastic case. To 


understand exactly why this is, we will proceed to fully carry out the derivation of 
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the stochastic SPE. The first derivatives in the new scales take the form 

d d(j) d dxo d dxi d 

dx dx d(j) dx dxQ dx dxi 

Id d d 
ed<i)^ dx.^^dx^ 
d dcj) d dxo d dxi d 

dt dt d(j) dt dxQ dt dxi 

- lA 

e (90 ’ 

where d(l)/dx = —1/e, dxo/dx = 1, dxi/dx = e, dxo/dt = dxi/dt = 0 and dcjy/dt = 
1/e. Note that we keep terms in the expansion up to the order of e. The second 
derivatives can now easily be written as 

d^ _ d^ d^ 2 (92 2 1 <9^ 

dx"^ dxi dxidxo e dxod(f> ^ dxi ‘^dxid(p e^ dcjP 
(92 _ 1 (92 
dU e^ (902 


( 6 . 20 ) 
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If the expansion (6.17) and the second derivatives (6.20) are inserted into the stochas¬ 


tic Maxwell equation (6.16), we obtain 


2 a^ 2 

a^o ^ ^^a^ia^o edxod(f)~^^ dxf ‘^dxid(f)~^ 

+ (U/. + .„) = Ig (at„ + + (6.21) 

1 a^ 

(a + i^^{x)) (eMo + Mi + ...) + 6—^^(eMo + Mi + ...)^ 

The terms up to 0{l/e) vanish. The terms of 0(1) yield 


a^^ a /aMo\ 
dxodcp dxQ \ d(p J 


( 6 . 22 ) 


This implies Mq is independent of xq, i.e., Mq = Mo(0, xi, X 2 ,...) and the derivative 
of Mq with respect to Xq goes to zero whenever it appears in the higher order terms. 
By taking this into consideration, we are left with the equation in the order of e 


a^Mi 

2_h 

dxiidcj) 


dxidcj) 


+ (a + vi^x)) Mo + 6 


d(fT 


(6.23) 


This is the hrst non-trivial order. The solution of this equation requires knowledge 
of Mq. If we eliminate the noise from this equation by setting z/ = 0, the solvability 
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condition of this equation leads to the deterministic short pulse equation (4.7), which 


we have derived in chapter four. 

The question then naturally arises as to how to obtain Mq in the presence of white 
noise. We hrst divide the evolution variable Xq into n equal periods such that the n 


periods are 0 to 1, 1 to 2, ..., Xn-i to Xn = xq. Let’s now integrate equation (6.23) 
with respect to Xq from zero to one; 


Iq dxod(j) Jo 


'd‘^Mo , 

2— -;:r“ + (ft + Z/^(x))Afo + b- 


dxid(f> 


d(fT 


dxn. 


(6.24) 


It is easy to differentiate the right-hand side because Mq is independent of Xq. Car¬ 
rying out the integration in both sides yields 


9Mi(0,xo = l,xi,...) 9Mi(0,xo = 0,xi, ...)- 


50 


50 


+ + i{x)dxo))Mo + h P [ dxo. 


d^Mo 


(6.25) 


5xi50 


502 


Note that changing the limits of integration from zero to another xq (other than one) 
does not make any difference in our discussion. We require the left-hand side to be 
zero because 5Mi/50 grows unbounded with time. In other words, the term in the 


left-hand side of equation (6.23) is a secular term, and as such it must be removed. 
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As a result, we pin down dMi/d(j) to zero at the end points of each period in order to 
avoid any growth of the function Mi. Notice that we can just set the left-hand side of 


equation (6.23) to zero without hrst integrating the equation. It must be emphasized 


that the integration of the noise term is to be handled carefully. When we integrate 


equation (6.23) from zero to one (over the hrst period), we draw a random number 


Cl = Jq ^(x) dxQ. In a similar manner, we integrate equation (6.23) from one to two 


(over the second period), kill the growth of Mi and draw another random number 
C 2 = Ji ^(x) dxQ. If this process is repeated n times, we obtain n random numbers 
such that (Cn) = (Ci)C 2,C35 ••■)Cn)- The collection of these random numbers produce 
a normal distribution as is the case with white noise. However, white noise is a 
continuous distribution, whereas the collection of these numbers is a discrete one. 
If we follow up on the discussion of the discrete approximation of continuous noise 


mentioned in the previous section, we obtain 


- 2(Mo)xi</, = (a uE{xi))Mo h{Ml)^^ . (6.26) 

This is the stochastic version of the short pulse equation corresponding to the deter¬ 


ministic form (4.7). If the transformation (4.21) [Aq is replaced with Mq) is applied to 
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equation (6.26), the stochastic equation corresponding to the deterministic equation 


of the form (4.8) can be obtained as 


UxT = il + ^SiT))U+^iU^)xx, 


(6.27) 


where we assume S(T) = S(—T). We address this form of the stochastic short pulse 
equation as the stochastic SPE. It is the governing equation for soliton propagation 
in the stochastic environment. 

It is noteworthy to mention that the strength of the slow-scale noise 5(a:i) is 
governed by the expansion parameter e as well as ly, whereas the strength of the fast 
noise ^(x) is controlled by z/. Since the leading order amplitude Mq is independent of 
Xo and the randomness is only dispersive, the strength of the slow noise rather takes a 
simple form. In cases where the leading order may be more complicated, the strength 
of coarse-graining noise (slow scale noise) requires a more careful treatment [53]. 


6.3 Numerical Analysis of the Stochastic SPE 

It is not a rhetorical question to ask how the ultra-short solitons of the deterministic 


short pulse equation propagate in a world that is not perfect. We have already 
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derived an equation to model the ultra-short pulse propagation in this imperfect 
world. The origin of the noise and randomness has already been mentioned in the 
previous section. However, we have not tested how these pulses evolve in a stochastic 
environment. In this section, we experiment with the pulse propagation in such an 
environment, and show some of our numerical results that analyze the evolution of 
the SPE solitons via the stochastic SPE. By doing so, we use the exact solution of 
the deterministic SPE as an initial condition and let it propagate in the stochastic 
short pulse equation. Figure |6T] shows the evolution of the SPE soliton at t = 51.2 
units. The dashed line is the solution of the stochastic SPE at t = 51.2, and the 
solid line is the analytical result at the same distance. The noise strength in this 
experiment is chosen as ly = 0.05. We set the soliton parameter m = 0.3 and the 
expansion coefficient e = 0.2 in a scheme that employs the semi-implicit method 
(to be discussed in greater detail in the succeding chapter). The noise is generated 
through the Matlab’s random generator rand in a normalized way. Although there 
is a small change in the shape of the soliton, it propagates stably in the stochastic 
SPE. The comparison of the solution of the stochastic equation and the exact result 
at t = 51.2 indicates that the effect of random dispersion on the ultra-short pulses is 
not strong . However, it is known that when the NLSE is used as a modeling equation 
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for very short pulses in optical fibers, the random dispersion is quite strong [5l] , and 
in that case the modihed NLSE must be used. On the other hand, the SPE solitons 
persist in the stochastic short pulse equation in the presence of noise. As long as 
the noise in the physical system does not embody large scale fluctuations, the SPE 
solitons undergo a stable propagation in a stochastic environment as conhrmed by 
our numerical results. 

In analogy with the deterministic case, we are interested in observing how the SPE 
solitons evolve in the stochastic nonlinear wave equation if they are used as initial 


conditions. Figure 6.2 displays the numerical result and its comparison to the exact 


result at the propagation distance t = 25.6 units. The dashed line is the numerical 
solution of the stochastic Maxwell equation, and the solid line is the exact solution 
without noise at t = 25.6. We use the Ablowitz-Ladik scheme in this experiment. The 
random numbers incorporated in our experiments are generated similar to the way in 
which we generate them for the stochastic SPE equation (refer to chapter seven). The 
noise strength z/, the soliton parameter m and the expansion parameter e are set to 
0.05, 0.3 and 0.2 respectively ffl- Note again that the initial condition and the exact 
result are modihed, respectively, as Uinuiai = eM(x/e,0) and Uanaiyticai = eu([(a: — 
25.6)/e], —e25.6) as imposed by the multiple scale expansion in the stochastic case 
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Figure 6.2: Evolution of the SPE soliton in the stochastic Maxwell equation. 
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like had been shown in the deterministic case. With the appropriate choice of initial 
condition, we observe in hgure [6^ a stable propagation of the solitons in the stochastic 
environment through the nonlinear wave equation as well. The comparison of the 
exact and numerical results clearly show that the impact of the noise somehow affects 
the soliton as it propagates, but these solitons persist in the stochastic nonlinear 
wave equation despite the presence of the of noise. On the other hand, it seems 
as though the impact of the noise on the solitons is more significant in the case of 
numerical solutions to the stochastic Maxwell equation than those for the stochastic 
SPE, which is used to model soliton propagation in the stochastic environment. In 
the next section, we will dicuss such a signihcance numerically and qualitatively. 


6.4 Comparison of the Stochastic SPE and Stochas¬ 
tic Maxwell Equation 


The stochastic SPE ( 6.27| ) (or alternatively (|6.26D) is expected to be a good approx¬ 


imation of the stochastic Maxwell equation (6.16) up to at least 0(l/e). We have 


already shown that the SPE solitary waves propagate in the deterministic nonlinear 
wave equation to a distance of 0(l/e^) ~ 26 (much larger than the leading order 
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distance) with very good agreement (see figure 4.4). The error accumulation versus 
propagation distance was also shown in figure |4.5 We now ask how much error is 
accumulated if we allow the SPE solitons to propagate in the stochastic Maxwell 
equation. Figure |6.3 demonstrates the error accumulations between the exact result 
(deterministic SPE solitons) and the results of the stochastic SPE and the Maxwell 
equation for one realization of the noise. For the simulation of the Maxwell equation, 
the error in the deterministic case is subtracted from the error in the stochastic case. 
The blue line shows the L°°-norm generated for the difference of the exact SPE soli¬ 
tons and the solution of the stochastic Maxwell equation with the deterministic error 


(L“-norm shown in figure 4.5) removed, whereas the green line is the L“-norm for 
the difference of the exact SPE solitons and the solution of the stochastic SPE. We 
observe an excellent agreement between the fast scale noise and the slow scale noise 
for this particular realization. 


Since we obtained figure 6.1 and figure 6.2 for one realization, repeating the ex¬ 


periments may result in a change in the shape of the pulses for other realizations due 
to the randomly distributed nature of the white noise in the system. Therefore, the 
error accumulation in a stochastic environment for one realization may be meaning¬ 


less. To compare the statistical properties of the fast noise and the slow noise, we 
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Figure 6.3: The growth of the deviations for the evolution of the SPE soliton in the 
stochastic SPE and stochastic Maxwell equation as dehned bv L°° norms. 
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have to draw many realizations and see whether or not the statistical properties are 
in good agreement. 


The path-wise correspondence in hgure |6.3| indicates that the related probability 
distributions must be in agreement with one another as well. In hgure |6^ we plot the 
probability distribution versus the L°°-norm of the difference between the stochastic 
PDEs and the exact SPE solitons at the propagation distance x = 25.6 units. The 
blue line is the probability distribution of the Maxwell equation and the green line is 
the probability distribution of the stochastic SPE. The plots are obtained for 10, 000 
realizations of the stochastic SPE and the stochastic Maxwell equation by joining the 
midpoints of the histograms of the deviations of the solutions to the stochastic PDEs 
from the deterministic evolution of the SPE solitons. The probability distributions 


match to very good accuracy as seen in figure 6.4, and this is indicative of a cor¬ 


relation between the noise over slow scales as well as fast scales in accordance with 
the discrete-noise approximation, coarse-graining noise and multiple scaling expan¬ 
sion (see chapter seven for more details). The tail of each distribution curve reflects 


the fact that any deviation between the numerical results and the exact result at 
X = 25.6, which happens to be greater than 0.02, becomes less probable [4U] . 
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Figure 6.4: Probability distribution of the deviations of the stochastic results from 
the deterministic results versus the L“-norm of the difference between the stochastic 
PDFs and the exact SPE solitons at the propagation distance x = 25.6 units 








Chapter 7 


Numerical Methods 


It is not rare to come across a problem whose analytical solntion may either not exist 
or be too tedious to obtain by hand in applied sciences and mathematics. Numerical 
analysis is the branch of mathematics and applied sciences that is used to find approx¬ 
imations to such difficult problems. These problems may include finding the roots 
of non-linear equations, solving differential equations, complex integration, numerical 
differentiation, Fourier analysis, hnite differences and so on [SSIE]. One can, on the 
other hand, use numerics as a tool to validate the analytical answers of mathematical 
objects and equations. In this sense, numerical methods are the experiments, and 
computers are the laboratories of mathematics. We are fortunate enough that there 
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are a variety of programing languages available to those which deal with such experi¬ 
ments in applied helds |56]. The software we have used extensively in the preparation 
of this thesis is Matlab. Matlab is a high-level technical computing language and 
interactive environment for algorithm development, data visualization, data analysis, 
and numeric computation |1]. One can either use Matlab’s built-in-commands to per¬ 
form tasks or write their own codes to execute [5]. Although compiled languages such 
as C, C-I--I-, and Fortran are the most efficient options in terms of execution speed, 
the most important advantage of Matlab is that it is a very simple, yet powerful, 
programming language. 

A number of different numerical schemes is widely available for solving numerical 
problems. In the effort of understanding the short pulse dynamics, we have used the 
schemes employing Euler’s Method, Runge-Kutta Method, the Semi-Implicit Method, 
Ablowiz-Ladik Scheme and the exponential time differencing method . Each scheme 
has its advantages and disadvantages in terms of numerical stability and computa¬ 
tional errors. We shall discuss, in this chapter, the numerical tools we have used in 
our codes and the application of each numerical method to the Maxwell and short 
pulse equations in both the deterministic and stochastic cases. 
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7.1 Fourier Transform and FFT Algorithm 

The Fourier transform is a powerful technique that applies to a wide variety of prob¬ 
lems in mathematics and applied sciences, and therefore, is widely used tool in nu¬ 
merical calculations as well. We discuss the details of the Fourier analysis, and how 
it is implemented in our numerical schemes in this section. 

The French mathematician, Fourier, found that any periodic waveform can be ex¬ 
pressed as a series of harmonically related sinusoids, i.e., sine and cosine waves, whose 
frequencies are multiples of its fundamental frequency or the hrst harmonic. The ba¬ 
sic idea behind this technique is to look at the problem from a different perspective. 
The periodic function is hrst transformed to a new space, called the Fourier space, in 
which it is represented as the sum of the sines and cosines. The manipulation of the 
function may be carried out with the new look of the function in the Fourier space, 
and a solution is sought. The Fourier transform is reversible and one can always 
transform the function back to the original space via the inverse Fourier transform 
once the mathematical operation(s) in the Fourier domain is done. 

Mathematically speaking, a periodic function f(x) can be expressed as a series of 
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the sines and cosines [571 EH] as 



n=oo 


n=oo 


( 7 . 1 ) 


i=0 


i=0 


where 



f{x)dx 


f {x)cos{nx)dx 


f{x)sin{nx)dx 


and n = 0,1, 2, 3,.... The hrst term ao/2 in the series is a constant and represents the 
average component of the function f{x). The terms with the coefficients a\ and hi 
in the series represent the fundamental frequency component w. Likewise, the terms 
with the coefficients 02 and 62 represent the second harmonic component 2 w, and 
so on. If the periodic function f[x) has the even symmetry, or in other words, if it 
is an even function, i.e., f{—x) = f{x), the series consists only of the cosine terms 
with zero or nonzero Qq. If it has the odd symmetry, that is, if it is an odd function 
(/(—x) = —/(x)), the series includes only the sine terms. 
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The Fourier series can be generalized to the complex numbers, and further gener¬ 
alized to derive the Fourier transform. We will not show the derivation of the Fourier 
transform here, but we will only give the dehnitions. The Fourier transform and the 
inverse Fourier transform are dehned respectively as 


F(w) 

fix) 


' —OO 

/*oo 


C^^f{x)dx 
e-^^^F{k)dk , 


(7.2) 


where w = 2nk is the angular frequency, k is the Fourier frequency and i = . 

The Fourier transform maps a time series into the series of frequencies of the am¬ 
plitudes and phases. The inverse Fourier transform maps the series of frequencies 
(their amplitudes and phases) back into the time series. The derivatives of the peri¬ 
odic functions can be also be written in the Fourier space. Let us show the Fourier 
transform of the hrst derivative as an example 


F'{w) = ffx) = / F^^f{x)dx 


— OO 
OO 


= f(x)e""‘= - iw / e"“’f(x)dx. 


( 7 . 3 ) 


— OO 


— OO 
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Since f(x) —)■ 0 as x —)■ ±oo, equation (7.3) becomes 


F'(w) = —iwF{w). 


We can generalize to the nth derivative; 


(7.4) 


F"(ui) = (7.5) 

Note also that the condition for the Fourier transform to be applied is that /(x) —?-0 
as X —)■ ±cx), and /(x) must be a periodic function. 

There is also the discrete counterpart of the Fourier transform [59l |57j and is called 
the discrete Fourier transform (DFT). The DFT can be turned into the numerical 
language easily and efficiently. Suppose we truncate the Fourier series at the Nth 
term and use the N number of harmonics; 

N-l 

= n = 0,l,...,iV-l (7.6) 

k=0 

This is the forward discrete fourier transform equation. The complex numbers /q, /i,..., /at 
are transformed into the complex numbers Fq, Fi,Fn- The backward formula, 
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which is called inverse discrete fourier transform (IDFT), can be written as 


iV 


A^-1 


^ — 2'Kikn/N 


n = 0,1, ■■■, N — 1 


(7.7) 


n=0 


The complex numbers Fq, Fi ,are transformed into complex numbers /q, /i ,/at 


by IDFT. These two formulae (7.6) and (7.7) are the basis of computer algorithms 


for the fourier analysis. 

The DFTs and IDFTs are computed using the so-called fast fourier transform 
(FFT) algorithms in modern numerical applications [60]. The FFT algorithms have 
been discovered independently by several researchers, and many people have since 
contributed to the development of the FFT schemes. The starting point for the 
modern usage of the FFT dates back to the seminal paper published by John Tukey 
of Princeton University and John Cooley of IBM Research in 1965. An FFT algorithm 
re-expresses the DFT of a size N = NiN 2 in terms of the smaller DFTs of sizes Ni and 
N 2 recursively to reduce the computation time. The FFT is a very efficient technique 
and replaces the DFT mainly because of two reasons such that the FFT generates very 


accurate results and is a quite fast algorithm. The computational time for the DFT 
algorithms is proportional to where N being the number of discretized points. On 
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the other hand, the FFT of a size N does Nlog 2 {N) number of operations to carry 
out the Fourier transform [59]. For example, if the number of data points is 1000, 
i.e., N = 1000, then the algorithm requires 1, 000, 000 operations to take the Fourier 
transform in case of using the DTF. However, an FFT based-algorithm would do 
approximately 10, 000 operations. This means that a hundred times less operations 
done by the FFT algorithm. This is an enormous computational cost saving. 

The derivation of the fast algorithm FFT starts with the dehnition of the discrete 
Fourier transform such that 


Fn 


N-1 

k=0 


N/2-1 


N/2-1 

+ /2fc+ie2™(2fc+l)/W 


fc =0 


k=0 


N/2-1 


I ^27vin/N 


N/2-1 


k=0 


k=0 


(7.8) 


This is the FFT algorithm. Notice that n runs from 0 to N, not just to N/2 in the 


last line of (7.8). The discrete Fourier transform of N length is reduced to the sum of 


two the Fourier transforms of lengths N/2. The reduction of each Fourier transform 


to a Fourier transform of a smaller size can be done recursively. Although there are 
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different adaptation of the FFT, the case where N equals to a power of 2 is especially 
attractive. If N is an integer number of power 2, the FFT of length N is expressed in 
terms of two FFTs of lengths N/2, then four FFTs of lengths N/4, then eight FFTs 
of lengths N/8 and so on until we obtain N numbers of FFTs of length one. An 
FFT of length one is just the number itself. If N is not an integer power of 2, it is 
still possible to express the FFT of length N in terms of the several shorter lenghts 
of FFTs. Breaking a big size of discrete fourier transform into a number of smaller 
sizes of DFTs has an enormous impact on computational time. For each value of n 


in (7.6), computation of requires N complex multiplications and — 1 complex 


additions. Therefore, computation of length N requires approximately N‘^ complex 
operations for big values of N whenever a DFT algorithm is utilized. However, if 


N = 2^ was chosen in (7.8), the number of steps in the recursion would be p. There 


are also N number of complex operations in the hnal stage of the FFT making the 
total number of computation Np = N{log 2 N) for the FFT algorithm. 

Let us mention how we implement FFTs and IFFTs in our numerical schemes 
before moving into the next section. We have extensively used Matlab’s fast Fourier 
transform algorithm by means of fft and ifft built-in commands in our research en¬ 
deavour. The built-in fft and ifft functions are based on the FFTW, the fastest Fourier 
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transform in the West, developed at MIT by Matteo Frigo and Steven G. Johnson. 
The periodic u fnnction of the SPE and Maxwell eqnation with vanishing bonndary 
conditions is manipnlated by forward and backward Fonrier transforms via fft and 
ifft commands. Each Fonrier transformed fnnction is mnltiplied by the angnlar fre- 


qnency. Since the short pnlse eqnation (4.8) and Maxwell eqnation (2.24) inclnde the 


Erst and the second derivatives of time, we also apply the derivative of the Fonrier 


transform to these eqnations according to (7.3) and (7.5) as well. The nnmerical 


schemes employing the nnmerical methods (except the ETD method) that will be 
discnssed in the next sections implement Matlab’s bnilt-in commands fft and ifft to 


solve the short pnlse eqnation and Maxwell eqnation. 


7.2 Solving Differential Equations: Euler’s Method 
and Midpoint Method 

The differential eqnations are commonly used for mathematical modelling of scientihc 
inquiries [2], and play a prominent role in many helds such as physics, engineering, 
and economics [T]. Whenever there is no analytical solution available, the numerical 
approximations are required. We will, in this section, discuss the Euler’s numerical 
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schemes, which we used them in the numerical analysis of Maxwell and the short 
pulse equations extensively. 

Let us mention Euler’s methods briefly Erst. Consider the hrst order differential 
equation 

I =/(..<). (7.9) 


Let [a, b] be the interval over which we want to hnd the solution to the above initial 
value problem (IVP) with a given initial value y{a) = y{0). Since we want to hnd 
a numerical solution, our aim is not to seek a differential function that satisfies the 
IVP. Instead, we will generate a set of points [{ynUn)] which will be used to satify 
the differential equation. Using the formal dehnition of the derivation and the Taylor 


expansion [55], one can approximate the left hand side of (7.9) for the small increment 
of time t, i.e.. At as 

Vn+l y-a 


dt 


t—tfi 


At 


where we use the notation yn for y (tn) and yn+i for j/(t„ + At). The differential 


equation (7.9) evaluated at time t = tn is then be 


Vn+l yn 

At 


/(t/n) ^n)) 


(7.11) 
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which can be rearranged to obtain Euler’s approximation 


Un+l Vn P f iVnytri) ■ 


(7.12) 


At t = 0, we have yi = yo + Atf{yo,to). Since the initial value ?/(0) is given at 
t = 0, one can calculate f{yo,to) and, therefore, the first iterated value yi. This 
iteration process is repeated until the iterated values approximate the solution curve 


y = y{t). This method is called the finite difference method, and equation (7.12) 


is the difference equation. For the reason that the difference equations approximate 
derivatives, the finite-difference methods approximate the solutions of the differential 
equations. 

The numerical techniques approximating the solutions of the differential equations 
may result in different results. Hence, error analysis is very important to see how 
good the numerical results are. Apart from having the round-off errors when using 
any finite difference methods to approximate the solution of the differential equations, 
there is also a discretization error and truncation error. The discretization error can 
be defined as the difference between the analytical solution and the numerical result 
obtained by the difference method. The truncation error arises from truncating the 
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Taylor series at the first derivative in case of Euler’s approximation (7.12). Both the 


discretization and truncation errors depend on the step size At or dt. The Euler’s 
method approximating the hrst derivative of a function accumulates truncation error 
in the order of the step size or discretized interval, i.e., 0(At). The smaller the 
discretized interval At is, the smaller the truncation and discretization errors are. The 


hrst order Euler scheme (7.12) may therefore give relatively big errors as the process 


proceeds, and therefore, it has limited usages. The error accumulation of the Euler’s 
scheme can also be qualihed if it is applied to the SPE. Forexample, when the initial 
ultra-short pulse obtained from equation (??) with m = 0.3 is allowed to propagate in 
the hrst order Euler scheme with a step size dt = 0.001, the result does not agree with 
the analytical solution with a good accuracy. The error accumulation is much bigger 
than the one obtained by Euler’s midpoint formula even at the shorter propagation 


distances as shown in hgure 7.1 The graph displays the diherence between the 


numerical and the analytical value of u (maximum error) at the propagation distance 


t = 15 (x is the temporal variable in the SPE (4.8)). 


We use Euler’s method in one step only to generate an initial value. Euler’s method 


(7.12) is implemented in the SPE numerics together with the midpoint method. The 
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Figure 7.1: The maximum error between the numerical solution of the Euler’s first 
order scheme and the analytical result (upper graph), and between the solution of 
the midpoint formula based scheme and the analytical result (bottom graph) at the 
propagation distance t = 15 units. 
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midpoint formula for numerical integration of (7.9) is 


2/n+l Vn—l A (jJnUn) ■ 


(7.13) 


The midpoint formula is the first central difference approximation and is a second 


order formula. It is just the improvement of the hrst order Euler formula (7.12) and 


can be derived in a similar manner (see the next section for the derivation). The 
midpoint formula may also be classihed as the second order Runge-Kutta formula 
ISH and is sometimes refered as the leapfrog method. The leapfrog method is widely 
used because of its good stability when solving the partial differential equations with 
the oscillatory solutions. The error at each step of the midpoint method is of the 


0{AP). As shown in hgure 7.1, midpoint formula generates much less error than 


Euler’s method does, and it may generate very good results at the expense of some 
more computational effort. 

The application of the Euler’s method and Midpoint method is straightforward 


in our numerical codes for the SPE. The Fourier transform of the SPE (4.8) with 


respect to x variable is hrst taken via fft command of Matlab, which then yields 
Ut = u/iw + [{iw)/C\u^. The zero frequency mode must be exluded in the above 
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equation because the first term in the right hand side otherwise becomes inhnity 
in the fourier domain. Once we hnish with the Fourier transform, we can take the 
inverse Fourier transform via ifft and be back in the spatial domain with only one 
t (evolution variable) derivative left, i.e., Ut = au + (3v?. The a and (3 are the new 
coefficents after the Fourier transforms being done. The iteration schemes can now be 
applied. The time and space domains can be chosen and discretized in a desired way 
in the codes. These choices have to be made by considering numerical stability and 


error accumulations. If n = 0 is chosen in (7.12), we have the hrst iteration equation 


ufei.ii) = m( 90, 4) + A«(ati(!/o,«o) + l3u(yo,to)'‘). 


(7.14) 


This is an initial value problem and one can use the analytical solution (4.27) at 


f = 0 for u{yoHo). Once the initial value is substituted into (7.14), the hrst iterated 


value is generated. This is the entire usage of the Euler’s hrst order formula in our 
numerical scheme. Following the application of the Euler’s method, we can now apply 


the midpoint formula. If we choose n = 1 in equation (7.13), we get 




(7.15) 
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Using the initial value and Euler’s result for u(yo,to) and respectively, we 

obtain the second value u{y 2 ,t 2 )- We get the third from the hrst and the second,the 
fourth from the second and the third, and so on by the successive application of 


equation (7.13). Note that Fourier tranform and inverse fourier transform is applied 


at each step. This iteration repeated until the process stops. Overall process of 
this scheme is a stable numerical propagation of the SPE solitions in the short pulse 
equation. Figures produced by midpoint formula are shown in chapter four already. 

At this point, one may argue the initial value obtained from the analytical solution. 
Note that the spatial variable t and the temporal variable x are the free variables of 


the partial differential equation (4.8). The initial condition (4.27) does not depend on 


the free variable x directly. Instead, it depends on another variable y from which we 
can get initial contion u{y, t). We cannot chose y values in an arbitrary way because it 


is not a free variable. It depends on x, t and itself, i.e., y{x, y, t) ( See equation (4.28) 


The hrst approach coming to mind might be writing y values in terms of x and 


t. A blind look at the equation (4.27) shows that it may not be possible to take the 


inverse of the equation (4.27). The y values are approximated from the free variables 


X and t numerically. To do this, we hrst set t = 0 and obtain y from equation (4.27) 
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as 


y = x- 


2mn{{msin{2ny) — nsinh{2my)) 
{rn?sin{nyY + rPcosh{myY)) 


(7.16) 


where m and n are the equation parameters and may take different values according 


to (4.28). We choose an interval for free x values and set y = x aX hrst. Having x 


and y values in hand by this choice will let us to apply equation (7.16) to get a new 


interval for y values. This new set of y values along with the x values can be used to 
generate another new set of y values. This process has been repeated until we generate 
a precise set of y values. How good the hnal values of y is the next question to ask. 


A simple way to check is to insert numerically produced y values in equation (4.27) 


and re-generate free x values from these y values. The plot y — x versus x in figure 


7.2 demonstrates that y and x values are not the same and, therefore, the necessity 
of the y values for the numerical experiments is not questionable. The logarithmic 


plot of the maximum error versus the number of iterations in hgure 7.2 displays that 


one has to repeat the iterations about sixty times to get precise values of y. Once the 
minimum number of iterations are executed, the numerical y values may be safely 
used in the experiments of the short pulse dynamics. In the case of approximate 


solution (4.30), such a problem does not arise because y ^ x. Before closing this 
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Figure 7.2: The plot of the difference between the numerical computation of the vari¬ 
able y and the free variable x versus x (upper graph), and the plot of the logarithmic 
difference versus the number of iterations (bottom graph). 
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section, we may make our final remark that our numerical scheme based on Euler’s 
formula and midpoint formula can be used for further research applications. 


7.3 Semi-Implicit Method and Euler Central For¬ 


mula 


The hrst numerical scheme adopted for the Maxwell equation (7.29) is the Euler’s 


central method. We have already discussed the hnite difference approximation for 
the hrst derivative in the previous section. In dealing with the SPE, we only need a 
difference formula based on the hrst derivative of the function u. However, the numer¬ 
ical technique needed for the nonlinear wave equation requires the second derivative 
of the function u{x,t). 

The derivation of the hnite diherence approximations for the derivatives of a 
function f{x) are based on forward and backward Taylor series expansions of f{x) 
about X [551 ESI [63]. We will only illustrate the derivation of the diherence formula 
for f {x) here. Let us start with Taylor expansion of f{x + h) 


f(x + A) = f(x) + h/(x) + ^f"(x) + ^f"\x) + .... 


( 7 . 17 ) 
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Similarly, 


f{x - h) = f(x) - hf’lx) + - ^/"(^) + 


(7.18) 


Here, h can be taken as the difference between the snccessive step in space or time 


domain, i.e., the step size. Adding np eqnations (7.17) and (7.18) eliminates the 


odd derivatives If we trnncate the series at the fonrth derivative, we 

obtain the relation for the second derivative f"{x) 


/"(x) = 


f{x + h)-2f{x) + f{x - h) 2h‘^f^^^{x) 




4! 


(7.19) 


The first term in (7.19) is the desired formula and called Euler’s central formula. The 


second term is the error term E{f, h) = {h‘^/12) due to the truncation. If we 


subtract equation (7.18) from (7.17) and truncate the series at the third derivative, we 


will obtain the midpoint formula (7.13) with the error term E{f, h) = (h^/6) / (x). 


Notice that the first order non-central forward finite difference (see equation (7.12)) 


can be obtained from equation (7.17) with the truncation error E{f, h) = {h/2) f”{x). 
The power of h in each error term shows the order of the truncation error for these 


hnite difference methods. As it is clearly seen in the formula (7.19), the formula 


is truncated at the fourth derivative and the error due to that is in the 0{h^). The 
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truncation error is not the only error source in the numerical work. We also encounter 
various kinds of errors when using an algorithm for the computations such as the 
discretization and round-off errors. The truncation error in computation arises from 
having a hnite number of terms instead of inhnitely many terms we have in theory. On 
the other hand, the round-off error is caused by storing numeric data in hnite bits, and 
the discretization error results from the fact that a function of a continuous variable 
is represented in the computer by a hnite number of evaluations. The discretization 
and truncation errors can usually be reduced by using a smaller step size. Such a 
decrease in the step size increases the computational cost in return. Note that the 
sum of the coefficients of / functions (f(x) and f{xEh)) is zero in all hnite diherence 
expressions. If h is very small, the values of f{x) and f{x Eh) will be approximately 
equal and the ehect on the roundoh error can be profound. On the other hand, we 
cannot make h too big because the truncation error would be intolerable. This means 
as we decrease the step size, we decrease the truncation and discretization error and 
increase the round-oh error. Therefore, we must hnd an optimum value of the step 
size which balances out these errors. One partial solution to this problem is to use 
a formula of higher order so that a larger value of h will generate a better accuracy. 
There may also be error due to the noise in the system and due to the instability of the 
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algorithm. Although we cannot get rid of these kinds of inevitable errors completely, 
we should minimize the error by choosing the optimum value of step size and by 
having a stable algorithm. The step size dt = 0.001 is in general used in the schemes 
employing the semi-implicit method and the Euler’s central formula. This choice of 
dt is small enough and computationaly affordable. 

Let us now look at the details of the scheme employing the Euler’s central for¬ 


mula (7.19) for the nonlinear wave equation (7.29). To apply the second order hnite 


approximation, we take the Fourier transform of equation (2.24) and apply the Euler 


central formula (7.19), which then leads to a difference equation in the Fourier domain 


as 


Un-\-l 2 'U 72 Ufi—\ -)“ Nt ( W T -)- ( W 


(7.20) 


where w is the angular frequency, n = 1,2,3...N and Un+i = u{wn,tn+i)- The choice 
n = 1 is the hrst iteration equation 


U2 = 2ui — Uq + 


{-w^ -F h)uo + (-C * 


(7.21) 


where U2 = u{w, 2At), ui = u{w, At) and uq = u{w, 0). One must have the knowledge 
of Ml and Uq to get the hrst iterated value U 2 - We use the Fourier transfom of the SPE 
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initial condition as the initial condition of the Maxwell equation in (7.21). This is done 


for two reasons; the analytical solution of the Maxwell equation is unknown as of yet, 
and the SPE is derived from the Maxwell equation. It may be worty to mention that 


the transformation rule (4.21) must be applied to the initial condition because this is 


not the exact solution of the nonlinear wave equation. If the SPE initial condition is 
not rescaled according to the transformation rule and is used as the initial condition 
of the Maxwell equation in the numerical experiments, the instablity is inevitable. 
The reader may ask the following question; how do we get hi? This is the other 
initial condition which includes the hrst derivative of u function. The easiest way of 
getting this second initial condition is to compute it numerically. The simple shift of 
h(tc,0) by At in the numerical code generates u{w,At). We also wanted to get hi 
analytically. This analytical implementation introduces a hrst derivative of u function 


(4.27). Since analytical solution (4.27) of the SPE does not depend on free variables 


X and t directly, a simple application of chain rule includes the time derivative of y 


function (7.16). After dealing with partial derivative of y, the hrst time derivative of 


u function can be written as 


— — -Ut + Uy{ — ft /(1 + fy)) 


( 7 . 22 ) 
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where Ut and Uy are the partial derivative of u with respect to t and y respec¬ 


tively, and ft and L are the partial derivatives of the second term in eqnation (4.27) 


{x = y + f{y,t)) with respect to t and y respectively. Once we get the hrst derivative, 
Ui can be obtained from the hrst order Euler’s method u {f) = [u{t + h) — u{t)]/h. 
Note that n = 0 choice instead of n = 1 as the starting value of n does not mean that 
we are stuck. In that case, backward Euler’s formula, u {t) = u{f) — u{t — h)/h, should 


be utilized to start the iteration scheme. The hrst iteration (7.21) with Fourier trans¬ 


form of Ml and Mo gives M 2 , which then together with mi generates M 3 . This iteration 
process is repeated many times until we hit the most precise value of the solution. 
Finally, we take the inverse Fourier transform of the latest iterated value and we have 
the numerical solution of u. As for the results produced by Euler central formula, hrst 


of all when the coefficients of nonlinearity in nonlinear wave equation (2.24) is taken 


0.005, it seems like the numerical scheme is doing a good job. Initial SPE solitary 
pulse propagates stably and is compared with the analytical solution very well at even 
long propagation distances such as 40 or 50 units distances. This way of comparing 
the numerical work with the analytical value can be misleading and deceptive. To be 
able to observe the ehect of nonlinearity with these coefficients of nonlinearity, one 
may need to propagate the initial pulse a much longer distances (forexample, 400 or 
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500 unit distances). Propagation to such a big distance will increase the computa¬ 
tional cost and computational time. Besides, numerical instability may interfere the 
numerical propagation for such a big propagation journey. Therefore, the coefficient 
of nonlinearity should be big enough in order to observe the effect of dispersion and 
nonlinearity within few units of propagation distance. Remember a soliton propa¬ 
gates stably because of the delicate balance between dispersion and nonlinearity. At 
the same time, we cannot choose these coefficients arbitrarily. For instance, if we 
just made both coefficients of dispersion and nonlinearity one, we would hurt such 
delicate balance between dispersion and nonlinearity, which then leads to an unstable 
propagation of solitary wave. For that reason, we rescale the coefficients of the non¬ 


linear wave equation (2.24) such that the coefficient of the dispersion becomes two. 


i.e., a = 2 and the coefficient of the nonlinearity becomes one third, i.e., b = 1/3. 
With these coefficients, we can now observe the effects of nonlinearity when we let 
the soliton progates a few units. As a matter of fact, this way of checking the propa¬ 
gation produces results within few units of evolution about the reliability of the code. 
The instability in the propagation of the initial pulse with these choice of coefficients 
starts even at two units of propagation slightly and builds up much more leading to 
the loss of the pulse as the pulse keeps moving more and more. It is not easy to 
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spot what causes this instability right away. Because there are two different ways of 
getting u(w, dt), and the usage of numerical shift u{w, dt) in the code (instead of the 
value u{w, dt) obtained by the analytical first derivative du/dt ) does not remove the 
instability in the code either, we cannot suspect wrong calculation of the analytical 


derivation of u function (4.27) with respect to t. Implementation of filtering was not a 


remedy to further application of our Euler’s central formula based numerical scheme. 
For this reason, another scheme is employed in the numerical analysis of the nonlinear 
wave equation. 

In the next section, we will discuss this new scheme, called Ablowitz-Ladik scheme. 
On the other hand, Euler’s central formula (midpoint formula) for the first derivative 
applies to the short pulse dynamics as discussed in the previous section and this 
implementation produces reliable results for spe. Before going into the next section, 
let’s also mention how it is implemented in a different way than the one mentioned in 
section 8.1. The central formula for the first derivative follows the basic manipulation 


of Taylor forward and backward manipulation. If the equation (7.18) is subtracted 


from the equation (7.17) and if the series is truncated at the third derivative, we get 




(7.23) 
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The dependence of trnncation error means that trnncation error here is in order of 


O(h^). The iteration equation follows the application of equation (7.23) to spe 


( 99 ) 


A simple treatment of a short pulse numerics using Euler’s central formula requires 
the fourier transform of the spe Erst. We can take care all integration in the fourier 
domain and transform back to the time domain after the integration is complete. For 
integration in Fourier space, we need initial contions u{yo,Wo) and u{yi,Wi) = Ui. 
Fourier transform of analytical solution at t = to gives u{yo,Wo) = Uq and application 


of (7.12) to Fourier tranform of the spe (4.8) gives u{yi, wi). Once we have the initial 


conditions, iteration scheme 


rUr. 


Un+1 = Un-1 + 2h[— + {iw)ul] 

XXL) 


(7.24) 


is ready to use. This formula is very similar to the equation (7.13). It may be called 


midpoint method in Fourier space alternatively. Here, n = 1,2,...,A^ and w is the 
Fourier frequency as mentioned before. Note that we can start n values from zero. 
In that case, we need to apply Euler’s backward formula to get m_i. When n = 1, we 


get an equation in Fourier domain similar to equation (7.15). Successive application 


of (7.24) results in «(?/„,Inverse Fourier transform of «(?/„, w„) would be the 
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final step in this scheme. The advantage of using this scheme is that it is much faster 
because Matlab does not apply Fourier and inverse Fourier transform n times, but 
only once. As for the results generated by these two schemes, they come out almost 
same. The ignorable difference between two results may stem from making coefficent 


of the fourier transform of spe (4.8) with respect to x zero whenever tc = 0 at each 


fourier step and/or some numerical uncertainty that we cannot state exactly. The 
other numerical technique that has extensively been used in the numerical analysis 
of the short pulse dynamics is semi-implicit method. It is the second order nonlinear 
method that is solving nonlinear partial differential equations as well as linear partial 


differential equations, and is simply an improved version of the Euler’s formula (7.12). 


There are also many other improved versions of Euler’s formula exist in the literature 
such as Runga-Kutta methods and Heun’s method [031133]. This modihed version is 


given as 


Un+l Vn _ r/Vri+l T Vn , \ 

At “ 2 ’ 


(7.25) 


where n = 0,1,2,...,A^ and At = h is the step size. The difference between semi- 


implicit method and Euler’s method can be seen by comparing the equations (7.25) 


and (7.12). When y and At are replaced with the short pulse variables u and h 
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respectively, and equation (7.25) is manipulated, the iteration scheme becomes 


_ I r p/'^n+l+'^n , \ 

^n+1 T hj ( — , j 


I rr /'^n+l+'^nx o/^n+1 + \3-| 

= Un + h[a{ ---) + /3{ ---) ]. 


(7.26) 


where a and (3 are the new coefficients as they appear in ut = au + (3u^. Note that 
t dependence in the second line in the above relation is embedded in the u function. 
The modihcation in the function f{u,t) (the hrst derivative of the u function with 
respect to t in our case) makes a scheme that works much better than the hrst order 
Euler’s scheme. A comparison between a partial modihcation and a full modihcation 
in the nonlinear short pulse equation demonstrates the degree of improvement in the 
numerical result. We obtain the scheme 


Un+l = Un + + ^ul]. 


(7.27) 


after we modify the linear part and leave the nonlinear part unchanged. Numerical 


schemes using iteration equations (7.26) and (7.27) solve the SPE (4.8) and the results 


are shown in hgure [Tffij An initial solitary wave of the short pulse equation propagates 












CHAPTER 7. NUMERICAL METHODS 


146 



X 


Figure 7.3: Comparison of the solutions of the semi-implicit formula to the exact 
result at f = 50 units of propagation distance. The solid line is the exact solution, and 
the dashed line displays the solution of the semi-implicit formula. The dotted line 
represents the solution of the numerical experiment with the semi-implicit method 
being only applied to the linear term of the short pulse equation. 
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50 units in these schemes and is compared with the analytical solution at t = 50 


units. Dashed line shows the result obtained by the scheme (7.26), whereas the dash- 


dotted line represents the result of (7.27) at t = 50. The solid line in hgure 7.3 


IS 


the analytical solution at the same propagation distance. Our numerical experiment 
indicates that the modihcation in both the linear and nonlinear parts approximates 
the exact solution almost perfectly. 

It may be useful to mention more about the details of the adopted semi-implicit 


scheme (7.26) in our numerics. The analytical SPE solution is used as the initial 


condition and the Fourier transform is applied to it before the semi-implicit iteration 


takes place. The choice of n = 0 in (7.26) leads to an 


. Ui+Uo. „,Ui+Uo-.3-. 

ui = uo + h[a {—-—) -F /f(—-—) J. 


(7.28) 


where uq is the initial condition and ui is the hrst iterated value. Unlike Euler’s 


hrst order method (7.12) and Euler’s central method (7.24), semi-implicit iteration 


equation (7.28) requires Ui as well as Uq as initial conditions to get the hrst iterated 


value. To resolve the problem of not having ui as the initial value, we choose ui = uq 


initially and apply semi-implicit method. The loop is executed three times to improve 
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the accuracy of the value Uj. After obtaining ui, the iteration is repeated for U 2 in the 


similar way. The process goes on until we obtain un- We choose the step size dt = 0.01 
in our code. The execution time is small. Although the step size is relatively large, 
the semi-implicit numerical implementation is fast and the highly accurate algorithm. 
Let us now compare the solutions of the semi-implicit formula with Euler’s central 


scheme (7.24). Figure 7.4 displays the difference of the numerical results obtained 


by these two schemes at t = 50. usemiimp is the numerical solution given by semi- 
implicit method and UEuierCentrai is the solution obtained by Euler’s central method. 


The numerical difference at f = 50 is almost ignorable as seen in hgure 7.4 It should 


be noted that both numerical results match with the analytical solution at f = 50 


almost perfectly (see hgure 7.3). Furthermore, the execution time of semi-implicit 
scheme to propagate the initial soliton hfty units is about four times faster than the 
execution time of Euler’s central scheme ( 38 seconds versus 135 seconds) to propagate 


it the same distance. 
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X 10"'’ 



X 


Figure 7.4: Maximum difference error between the numerical solutions obtained by the 
semi-implicit formula and the Euler’s central formula at t = 50 units of propagation 
distance. 
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7.4 Ablowitz-Ladik Method 

The Ablowitz-Ladix (AL) scheme is a fast and efficient nnmerical method for nonlinear 
evolntion eqnations snch as the nonlinear Schrodinger, Korteweg-deVries and modified 
Korteweg-deVries eqnations. For that reason we apply this scheme for the nonlinear 
Maxwell eqnation. 

Let ns start onr discnssion by hrst introdncing the AL scheme. Consider the 


nonlinear maxwell eqnation (2.24) in the form 


utt - u^x = F{u{xU )), 


(7.29) 


where F{u) = hu+c{v?)tt , and b and c are constant coefficients. To solve this eqnation 
nnmerically, we approximate all the derivatives by hnite differences. The domains in 
space and time are partitioned nniformly so that we have a mesh xo,xi, ...,xj for 
space and a mesh to, fi,...., fw for time. The difference between two consecntive space 
points is dx and between two consecntive time points is dt. The solntion u{x,t) at 
space point xj and at time tn along with space and time discretization can be written 
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as 

Xj = xo + dx, j = (0,1, 2...J) 

tn = to + dt, n = (0,1,2...A^) (7.30) 

u] = u(Xj,tn) 

We write finite difference second order derivatives (central differece formulas) for the 
space at position Xj and for the time at time respectively as 


^ - 2u] + 

dx"^ dx"^ 

d^u _ - 2u^ + 

dC dU 


(7.31) 

(7.32) 


Once we substitute equations (7.31) and (7.32) into equation (7.29) and rearrange 


terms, we obtain the following equation for 




+ + u]_,) + 2(1 - U)u] + dHF{u]) 


(7.33) 


This equation can be stabilized by setting r = 1 and using an average of the spatial 
coordinates for the function F{u) as shown by Ablowitz, Kruskal and Ladik [M] and 
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such a choice leads to the hnal form of the equation 


= -«r‘ 


+ T "Wj-i + dt^F[ 


u]+i+ul-i 


(7.34) 


This is the Ablowitz-Ladik (AL) scheme. This scheme has been tested on traveling 
wave and periodic breather problems over long time intervals and has given good 
results in terms of computational accuracy and computational costs |65l [66]. The 
comparison between the Ablowitz-Ladik scheme and the hnite difference schemes 
for the nonlinear Schrodinger equation indicates that the Ablowitz-Ladik scheme is 
faster than the latter [66|. We adapt the AL scheme for the Maxwell equation. The 
implementation of the AL scheme requires an initial condition. The SPE soliton 
solution is used as the initial contion. As discussed earlier, the SPE soliton stably 
propagates in the nonlinear wave equation, and therefore, this choice is not random. 
Nonetheless, we have to modify this initial condition acording to the multiple-scale 
parameters. The parameter e, which arises in the derivation of the SPE from the 


wave equation (see (4.1)), is chosen 0.2. As it is forced by multiple-scale expansion 
technique, e must be a small number and much less than one. There is no specihc 
reason why we choose e = 0.2 other than e being just a small number. On the 
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otherhand, one may choose e = 0.1 and should generate the similar results we got 


with e = 0.2. We also have to modify the spatial and time domains according to (4.2). 


The numerical scheme without proper scaling makes initial spe solitary wave condition 
spreaded around as it propagates. With the right way of choosing initial conditions, 
the Ablowitz-Ladik numerical scheme is stable. As a stability requirement, we choose 
dx = dt and equate them to 0.0125 in our code. The step size in this scheme is 
almost ten times bigger than the choice of the step size done in the Euler schemes. 
Such a change in the step size would effect the computational time incredibly if the 
Euler’s schemes are used. However, changing the step size from 0.0125 to 0.0031 
(approximately four times smaller now) in the AL code does not really alter the 
result, but increase the computational time only about twenty times more. Most 
importantly, the initial solitary wave propagates stably in the AL scheme unlike in 
Euler’s central scheme. The AL schemes’ results compare with the analytical solution 
at a very good accuracy at any propagation distances. We must note that the hltering 
in the Fourier domain removes the noise in the numerical scheme. The scheme is not 
stable without the hltering. Filter in the Fourier domain applies at every step of the 
iteration. We iterate once, and then we hlter out. This is the way we go until the 


very end. 
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The first iteration equation follows the choice of n = 0 and j = 0; 


4 = 

= —Uq^ + 1 + 


■U? + 1 


+ P 


■U? + ^ 


, (7.35) 


where a and /3 are constant coefficients of the Fourier tranformed SPE. The terms 
and u^i are the initial conditions obtained from the exact solution such that 
= u{x,t = 0) and = u{x, —dt). As for the term it is obtained by shifting 
the exact solution by one temporal step, i.e., Mq ^ = u{x — dx, t). In a similar fashion, 
the next iteration formulae can be obtained and used to hnish the integration. 


7.5 Exponential Time Differencing Method 

The exponential time differencing (ETD) numerical technique is another powerful 
method that solves nonlinear partial differential equations. We implement the modi- 
hed exponential time-differencing fourth-order RungeKutta method to solve the short 
pulse equation. Although the schemes employing the semi-implicit and Euler central 


methods produce very reliable results for the short pulse dynamics we have discussed 
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so far, they are not stable enough to test the particle properties of the SPE solitons, 

i.e., collisions of the SPE solitons. We implement the ETD method to allow the SPE 

solitons collide. The results presented in chapter Eve regarding the collision of the 
two -soliton solution of the SPE employ the ETD method. We will only discuss the 
ETD method itself in the rest of this section. 

The short pulse equation can be given in the general form 

■Ut = Lm + N(m, f), (7.36) 

where L and N are the spatially discretized linear and nonlinear operators respec¬ 
tively. Let us define 

V = e-^^u . (7.37) 

The factor is called the integrating factor. The time derivative of the integrating 
factor gives 

Vi = e-^*N(e-^'t;). (7.38) 


For the time stepping of the v function, the fourth order Runge-Kutta method is 
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used and given as 


a = dtf{Vn,tn) 


b = dtf{vn + a/2, tn + dt/2) 

c = dtf{vn + b/2,tn + dt/2) (7.39) 

d = dtf{Vn + cHn + dt) 


Vn +1 — VnP — (a + 26 + 2c + d), 
6 


where 


(7.38). 


dt is the time-step (discretization) and / is the right-hand side of equation 
If the V function is integrated over a single time step dt, we get 


f*dt 


Un-\-l — 6 Ur, 


+ e 


"Ldt 


e N{u{tn + T),tn + T)dT . 


(7.40) 


The integration on the right-hand side of the above equation has to be discretized so 
that the proposed generating formula [67] becomes 

s—1 m /■ \ 

Un+l = e^‘^^Un + dm y^(-l)^ f ^ j N„-fc , (7.41) 

where s is the order of the scheme. When the integration is carried out, the complex 
analysis is used to compute the coefficients via contour integrals in the complex plane. 
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We use the fourth order Runge-Kutta method so that s = 4 in our case. Therefore, 
the numerical iteration formula becomes 


ar^ = + L-i n(u„, U) 

+ L-i ^ ^^/2) 

+ L-i ^2N(6„, + dt/2) - N(u„, f,)) 

[-4 - Ldt + (4 - 3Ldf + (Ldtf)] N(m„, 

+ 2 [2 + Ldf + (—2 + Ldf)] (N(a„, + df/2) + N(6„ + + (if/2)) 

+ [-4 - 3Ldf - (Ldtf + (4 - Ldt)] N(c„, + dt)^ . 

(7.42) 


This is the iteration formula we use in our numerics. The initial condition is chosen as 


the two-soliton solution obtained from equations (4.35), (4.36) and (4.37) by setting 


W = 4 to produce figures 5T and 5T in chapter hve. The number of points on 
the contour integral and the time-step were chosen, respectively, as m = 32 and 
dt = 0.0125 in producing these hgures. The maximum error (4.913710“®) was also 
shown in chapter hve for the choices of the number of discretization N = 2^® and 


the time-step dt = 0.0125. The maximum error, on the other hand, are 0.0020939 
and 0.027373 units for the choices of iV = 2^^ and dt = 0.025, and N = 2^® and 
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dt = 0.05 respectively. The execution time for the choice of iV = 2^® and dt = 0.0125 
is about 1568 seconds, whereas the execution time for N = 2^^ and dt = 0.025 
becomes only 330 seconds. The numerical experiment in hgure [TTS] is the repetition of 



Figure 7.5: Interaction of the two-soliton solution of the SPE with the time step 
dt = 0.025 and the number of discretized temporal interval N = 2^^ 


the experiment of hgure 5T with the different discretized temporal domain and the 
time step, i.e., N = 2^^^ and dt = 0.025. Although the maximum error is bigger in 


this case, the results seems to be approximating the true values at a good accuracy 
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as seen in figure 7.5 If we set N = 2^^ and dt 


0.05, the integration time decreases 



Figure 7.6: Interaction of the two-soliton solution of the SPE with the time step 
dt = 0.05 and the number of discretized temporal interval N = 2^^ 


to only 82 seconds. However, hgure 7.6 indicates that the scheme is not stable. We 


observe some noise in the right tail of the slower soliton after the interaction. The 


noise like fluctuations shall be associated to the numerical instability. . 
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7.6 Random Numbers and Random Number Gen¬ 
erators 

The schemes approximating the solutions of the stochastic short pulse equation and 
stochastic nonlinear wave equation implement random numbers to account for the 
random variable and randomness in the equations. This section describes random 
number generators and how randomness is implemented in our numerical methods 

[6311591168]. 

A random variable is a variable whose possible values are numerical outcomes of 
a random phenomenon. There are two types of random variables; discrete random 
variable and continuous random variable. Numerical methods utilize the discrete one. 
A discrete random variable is a random variable that is countable and takes discrete 
values such as 0,1,2,3,.... We use simulated random variables to explain statistical 
pattern recognition and measures. The ability to generate random variables from 
known probability distributions is the subject of the computational statistics [69] . 
The statistical analysis of a random process in which there is a random variable that 
exhibits stochasticity and probabilistic features is called a probability measure. We 
have already shown probability measures and statistical features of the stochastic 
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short pulse (6.26) and Maxwell (6.16) equations in chapter six. When dealing with a 


random dataset related to the evolution of a random variable, we consider such dataset 
as one realization of an ensemble that consists of a large number of realizations of 
a generating process. This is the so-called random process or stochastic process. A 
random dataset displays the values of random variables. Such random values may 
be generated by several different ways. In our numerical methods, we use Matlab’s 
built-in random number generator randn that generates random numbers having a 
normal (Gaussian) distribution with a mean of zero and a variance of unity. In other 
words, randn returns a scalar value drawn from a normal distribution with mean zero 
and standard deviation one. The randn generator should not be mixed up with the 
rand generator which produce random numbers having uniform distribution. Note 
that the command rand returns a random number between zero and one. The choice 
between the usage of these two generators is related to the physics of a problem. 
As mentioned in chapter six, we model the noise in the system with the white noise 
whose statistical properties obey normal distribution. The Gaussian white noise is 
a good approximation of many real-world situations and provides maneuverings for 


mathematical models. 




CHAPTER 7. NUMERICAL METHODS 


162 


The numbers generated by Matlab’s generators rand and randn are not truly ran¬ 
dom numbers. They are, instead, pseudorandom numbers. Pseudorandom numbers 
generated by the random number generator is a sequence of numbers that approxi¬ 
mates the properties of random numbers. This means there is a period of repeating 
the sequence of random numbers. The sequence of numbers produced by Matlab’s 
pseudorandom generator randn is determined by the internal state of the generator. 
There are two methods, state and seed, that determines the internal state. The 
method state uses Marsaglia’s ziggurat algorithm, which is the default in Matlab’s 
versions 5 and later. The period for state generator is approximately 2®“^. The other 
method seed uses the polar algorithm, which is the default in Matlab’s version 4 and 
its period of repeating the pseudorandom numbers is approximately (2^^ — 1) * (vr/S). 
One may set the internal state to either one. However, changing states does not im¬ 
prove any statistical properties. Furhermore, randn will generate the same sequence 
of numbers in each session unless the state is changed since Matlab resets the state 
at start-up [1]. 

When using the randn generator with the either state in the stochastic SPE 
solvers, one must be very cautious. First of all, the period of the states may seem 
good enough, but it may not be. The generator randn may draw big random numbers 
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quickly from the normally distributed pseudonumbers. The impact of these relatively 
big numbers may diverge the pulse propagation since these numbers would be in¬ 
terpreted as a strong noise pertubations. Secondly, the discrete noise approximation 
necessary to implement the numerical method will modify the variance of the randn 
generator and, thirdly, the strenght of the noise will be modihed as imposed by the 
multiple scale expansion. Therefore, we have to multiply the randn generator by the 
coefficient eyjNoiseAmplitude.dt. The coefficient will surely take different values for 
the different amplitudes of the noise strength (noise amplitude) and dt. One may 
choose the noise amplitude and dt in a way that the relatively big random numbers 
drawn by the randn are yet to be rescaled properly. 

There are many different way of generating pseudorandom numbers [69]. We also 
implement rand generator in our numerical schemes. We force the rand generator to 
produce negative and positive pseudorandom numbers between zero and one. Such an 
implementation serves as a safety net to block the any unwanted big pseudorandom 
number in the system’s evolution. It may be worth to mention that the mean of 
these pseudorandom numbers may not most likely be zero. We produced the hgures 
in chapter six regarding the stochastic dynamics using the forced rand generator. 
The results clearly show the impact of the stochasticity on the ultra-short solitons. 



Chapter 8 


Conclusion 


The short pulse equation (SPE) possessing the exact solitary wave solution may be a 
better modelling equation for the ultra-short pulse propagation in the optical hbers. 
Our numerical work validates the exact derivation of the analytical solution of the 
SPE. We have also numerically shown that the SPE solitons approximates the solution 
of the Maxwell equation in one dimension. The higher order SPE is derived. The 
numerical validation of the higher order SPE and its comparison to the nonlinear 
wave equation are yet to be done. 

The solitonic properties of the SPE solitons have been tested by allowing them to 
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collide numerically. The particle-like behaviour of the SPE solitons shall be experi¬ 
mented in the presence of stochastic perturbations. 

The SPE is a one dimensional model. The short pulse equation in at least two 
dimensions may be more appealing. Although it seems that the multiple scale ex¬ 
pansion ansatz used in the derivation of the one dimensional SPE may not be the 
multiple scales for a two or three dimensional modelling equation, there may still 
exist the two or three dimensional SPE. 

Stochastic perturbations are present in all physical systems. The randomness has 
been introduced in the dispersion coefficient of the Maxwell equation based on the 
arguments made earlier. We have derived a stochastic version of the SPE and shown 
the results of the numerical experiments to underpin the impact of the stochastic 
perturbations on the SPE solitons. The future pursuit of the short pulse dynamics in 
the presence of the randomness may include the introduction of the stochasticity in 
the nonlinearity. In this case, the study of the statistics of the coarse-graining noise 
will require a more cautious approach. 

The transformation rule between the SPE and the sine-Gordon (sG) equation 
opens a channel by which the exact kink and anti-kink solutions of the sG equation 


are used to obtain the exact solutions of the SPE. The same transformation rule 
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may be applied to the stochastic sG equation to derive a stochastic form of the SPE. 
Whether the randomness is introduced in the dispersion or in the nonlinearity of the 
Maxwell equation, the stochastic SPE and its counterparts, if they so happen to exist, 
may be compared to the other stochastic models such as the stochastic NLSE and the 
stochastic KdV. The higher order stochastic SPE may be derived and experimented 
upon so as to improve the numerical results pertaining to the leading order stochastic 
SPE. 
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